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ABSTRACT 

We present the first secondary eclipse and phase curve observations for the highly eccentric hot 
Jupiter HAT-P-2b in the 3.6, 4.5, 5.8, and 8.0 /im bands of the Spitzer Space Telescope. The 3.6 
and 4.5 /im data sets span an entire orbital period of HAT-P-2b (P = 5.6334729 d), making them 
the longest continuous phase curve observations obtained to date and the first full-orbit observations 
of a planet with an eccentricity exceeding 0.2. We present an improved non-parametric method for 
removing the intrapixel sensitivity variations in Spitzer data at 3.6 and 4.5 fim that robustly maps 
position-dependent flux variations. We find that the peak in planetary flux occurs at 4.39±0.28, 
5.84±0.39, and 4.68±0.37 hours after periapse passage with corresponding maxima in the planet/star 
flux ratio of 0.1138%±0.0089%, 0.1162%±0.0080%, and 0.1888%±0.0072% in the 3.6, 4.5, and 8.0 ^^m 
bands respectively. Our measured secondary eclipse depths of 0.0996%±0.0072%, 0.1031%±0.0061%, 

0.071%+°;°^^^^, and 0.1392% ± 0.0095% in the 3.6, 4.5, 5.8, and 8.0 fiin bands respectively indicate 
that the planet cools significantly from its peak temperature before we measure the dayside flux 
during secondary eclipse. We compare our measured secondary eclipse depths to the predictions 
from a one-dimensional radiative transfer model, which suggests the possible presence of a transient 
day side inversion in HAT-P-2b's atmosphere near periapse. We also derive improved estimates 
for the system parameters, including its mass, radius, and orbital ephemeris. Our simultaneous fit 
to the transit, secondary eclipse, and radial velocity data allows us to determine the eccentricity 
(e = 0.50910 ± 0.00048) and argument of periapse {uj = 188.09° ± 0.39°) of HAT-P-2b's orbit with 
a greater precision than has been achieved for any other eccentric extrasolar planet. We also find 
evidence for a long-term linear trend in the radial velocity data. This trend suggests the presence of 
another substellar companion in the HAT-P-2 system, which could have caused HAT-P-2b to migrate 
inward to its present-day orbit via the Kozai mechanism. 

Subject headings: planets and satellites: general, planets and satellites: individual: HAT-P-2b, tech- 
niques: photometric, methods: numerical, atmospheric effects 
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1. INTRODUCTION 



The supermassive [Mp — 9 Mj) Jupiter sized {Rp 
1 Rj) planet HAT-P-2b (aka HD' 147506 b) was first dis- 
covered from transit observations by 



using the HATNet (Bakos et al. 20' 



^Bakosern!l(|2007b| 
02[ |2004[ ) net^;F5?ra 



ground-based telescopes, l^bilow-up radial velocity mea- 
surements of the HAT-P-2 system revealed that the or- 
bit of HAT-P-2b is highly eccentric (e = 0.5 Bakos et al 
2007b[ [Loeillet et al.||2008[ ). Only a handful ot transitmg' 



exoplanets have been shown to posses eccentricities in ex- 
cess of that of Mercury (e = 0.2), which has made HAT- 
P-2b an interesting target for many theoretical studies 
concerning the evolution of the HAT-P-2 system ([Jack- 
son et a l. 20 08; Fabryck;y[|2008 l [Matsumura et al.||i^008l 
Baraffe et""arT |2008l) . Because of its mass, HAT-P-2b 
lass or 



represents a class ot exoplanets that provides an impor 
tant link betw een extrasolar giant planets and low-mass 
brown dwarfs (Baraffe et al. 2008). Observational con- 
straints on the basic atmospheric properties of HAT-P-2b 
will provide an important probe into the structure and 
evolution of not only HAT-P-2b, but an entire class of 
massive extrasolar planets. 

Atmospheric circulation models for planets on eccen- 
tric orbits show significant variations in atmospheric 
temperature and wind speeds that provide an impor- 
tant probe into atmospheric radiative and dynamical 
timescales (Langton & Laughlin '2008; Lewis et al.[[2010 



Cowan fc Agol 2011- Kataria et ai.^ 20ll| F 
dent Mux on HA'i'-P-2b trom its stellar host 



The mci- 
host at periapse 
is ten times that at apoapse, which should cause large 
variations in atmospheric temperature, wind speeds, and 
chemistry. Heating and cooling rates for HAT-P-2b can 
be constrained by measuring planetary brightness as a 
function of time. Such observations are analogous to pre- 
vious phase curve observations of HD 189733b using the 
Spitzer Space Telescope, which provided the first clear 
observational evidence for atmospheric ci rculation in an 
exoplanet atmosphere (Knutson et al. 2007^ 2009b^|2012 ). 

HAT-P-2b's large orbital eccentricity makes it an ideal 
target for investigating hot Jupiter migration mecha- 
nisms. Gas giant planets such as HAT-P-2b are expected 
to form beyond the ice line in their protoplanetary disk 
far from their stellar hosts. HAT-P-2b currently resides 
at a semi-major axis of 0.07 AU from its host star, indi- 
cating that it must have migrated inw ard via a process 
such as gas disk migration (Lin et al. 1996j ), planet-pl anet 
scattering ("Ra sio fc Ford|l55t)[ ), secular i nteractions ([Wu 
fc Lithwick 201ll ), or Kozai migr ation ( Weidenschilling 
^Marzari||19g5nNaoz et al.[[2011[ ). HAT-P-2 b's close-in 
and highly eccentric orbit tavors one of the latter three 
mechanisms since disk migration tends to damp out or- 
bital eccentricities. For Kozai migration, the presence of 
a third body with at least as much mass as HAT-P-2b 
is needed. In this study we present six years of radial 
velocity measurements for this system, which allow us to 
search for the presence of a massive third body in the 
HAT-P-2 system. 

Here we present our analysis of the 3.6 and 4.5 /im 
full-orbit phase curves of the HAT-P-2 system, which 
include two secondary eclipses and one transit at each 
wavelength. These full-orbit phase curves represent the 
longest continuous phase observations ever obtained by 
the Spitzer Space Telescope. The orbital period of HAT- 



P-2b (5.6334729 d) is more than 2.5 times that of other 
exoplanets with pu blished full-orbit ph ase curve observa- 
tions: WASP-12b (jCo wa n et al.|2012a[ ) and HD 189733b 
(Knutson et al.[[20l^). Additionally, we present an anal- 



ysis ot previous partial orbit phase curve and secondary 
eclipse Spitzer observations at 8.0 and 5.8 fim respec- 
tively. We use these observations to characterize the 
changes in the planet's emission spectrum as a function 
of orbital phase and to probe the atmospheric chemistry 
and circulation regime of HAT-P-2b. The following sec- 
tions describe our observations and data reduction meth- 
ods (Sj2| and the results of our analysis ([js]). Section |4] 
discusses the results from our analysis of the Spitzer 
data and compares them to predictions from atmospheric 
models for HAT-P-2b. Additionally, we discuss trends in 
our radial velocity data that indicate the presence of an 
additional body in HAT-P-2 system in Section [4| Sec- 
tion [5[ overviews the main conclusions from our analysis 
and presents ideas for future work. 

2. OBSERVATIONS 

We analyze nearly 300 hours of observation of HAT- 
P-2 at 3.6 /im and 4.5 /im obtained during t he post- 
cryogenic m ission of the Spitzer Space Telescope (Werner 



et al.|2004[ ) using the IRAC instrument ( [Fazio et a"i.|2004l) 
in subarray mode. The observation periods were U'l' 
2010 March 28 to UT 2010 April 03 and UT 2011 July 
09 to UT 2011 July 15 for the 3.6 and 4.5 /im bandpasses 
respectively. Both observations cover a period just over 
149 hours with two approximately 2-hour breaks for data 
downlink, corresponding to ^1.2 million images in each 
bandpass. Each observation begins a few hours before 
the secondary eclipse of the planet, continues through 
planetary transit, and ends a few hours after the subse- 
quent planetary secondary eclipse. 

We also analyze observations of the HAT-P-2 system at 
5.8 iim and 8.0 /im obtained during the cryogenic phase 
of the Spitzer Space Telescope mission. The 5.8 /im ob- 
servations cover a single secondary eclipse of the planet 
that occurs on UT 2009 March 16. The 8.0 /im observa- 
tions cover a portion of the planet's orbit that includes 
transit, periapse passage, and secondary eclipse on UT 
2007 September 10-11. The 5.8 /im observations were 
obtained using subarray mode, while the 8.0 /im obser- 
vations were obtained using the full IRAC array. We 
obtain our 3.6, 4.5, 5.8, and 8.0 /im photometry from 
the Basic Calibrated (BCD) files produced by version 
18.18.0 of the Spitzer analysis pipeline. 

In subarray mode, 32x32 pixel images are stored in 
sets of 64 as a single FITS file with a single header. 
We calculate the BJD_UTC at mid-exposure for each 
image from time stamp stored in the MBJD_OBS key- 
word of the FITS header, which corresponds to the start 
of the first image in each set of 64. We assume uni- 
form spacing of the images over the time period defined 
by the AINTBEG and ATIMMEEND header keywords. 
The image spacing is roughly equal to the 0.4 s expo- 
sure time selected for the 3.6, 4.5, and 5.8 /im obser- 
vations of HAT-P-2 {K,r,ag = 7.60). For the 8.0 /im 
full array observations, the MBJD.OBS, AINTBEG, and 
ATIMMEEND header keywords are used to calculate the 
BJD_UTC at mid-exposure for each image, which are 
spaced by roughly the 12.0 s exposure time. In order to 
convert from UTC to TT timing standards as suggested 
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by Eastman et al. (2010), 66.184 s would be added to 
our BJD_U'1'C values. We report all of our timing mea- 
surements using BJD_UTC for consistency with other 
studies. 

2.1. 3.6 and 4-5 ^im Photometry 

We determine the background level in each image from 
the region outside of a ten pixel radius from the cen- 
tral pixel. This minimizes contributions from the wings 
of the star's point spread function while still retaining 
a substantial statistical sample. We iteratively trim 3ct 
outliers from the background pixels then fit a Gaussian 
to a histogram of the remaining pixel values to estimate 
a sky value. The background flux is 0.6% and 0.2% of 
the total flux in the science aperture for 3.6 and 4.5 /^m 
observations respectively. We correct for transient hot 
pixels by flagging pixels more than 4.5(t away from the 
median flux at a given pixel position across each set of 
64 images then replacing flagged pixels by their corre- 
sponding position median value. 

Several different methods exist to determine the stellar 
centroid on t he Spitzer arra.y suc h as flux-weighted cen- 
troiding (e.g. Knutson et al.|200"8 ), Gaussian centroiding 



and least asymme try methods (e.g. Stevenson et al. 2010 
Agol et al. 2010). We compare photometry calculatec 



using both Gaussian and flux- weighted centroiding esti- 
mates and find that for the HAT-P-2 data flux-weighted 
centroiding produces the smallest scatter in the final light 
curve solutions. This is in contrast with the w ork by 
Stevenson et all (120 101) and I Agol et all pOTol), which 



advocate Gaussian fits and least asymmetry methods to 
determine the stellar centroid for observations lasting less 
than 10 hours. We find that flux-weighted centroids give 
more stable position estimates over long periods, espe- 
cially if the stellar centroid crosses a pixel boundary dur- 
ing the duration of the observation. For these data, we 
calculate the flux-weighted centroid for each background 
subtracted image using a range of aperture sizes from 2.0 
to 5.0 pixels in 0.5 pixel increments. We flnd that the 
stellar centroid aperture sizes that best reduce the scat- 
ter in the final time series are 4.5 and 3.5 pixels for the 
3.6 /im and 4.5 fim observations respectively. 

We estimate the stellar fiux from each background sub- 
tracted image using circular aperture photometry. A 
range of aperture sizes from 2.0 to 5.0 pixels in 0.25 pixel 
increments were tested to determine the optimal aper- 
ture size for each observation. Additionally, wc tested 
time-varying aperture sizes based on the noise pixel pa- 
rameter described in the Appendix We find that we ob- 
tain the lowest standard deviation of the residuals from 
our best-fit solution at 3.6 /im using a variable aperture 
size given by the square-root of the measured noise pixel 
value for each image (median aperture size of 2.4 pix- 
els). For the 4.5 /im observations we find that a fixed 
2.25 pixel aperture gives the lowest scatter in the flnal 
solution. We remove outliers from our flnal photometric 
data sets by discarding points more than 4.5ct away from 
a moving boxcar median 16 points wide. We find that 
using a narrower median filter or larger value for the out- 
lier cutoff will miss some of the significant outliers. We 
also find that using a wider median filter or smaller value 
for the outlier cutoff will often selectively trim the points 
at the top and bottom of the 'saw-tooth' pattern seen in 
the resulting photometry for the 3.6 and 4.5 /xm observa- 
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Figure 1. Measured x (a), y (b), noise pixels (c), and raw pho- 
tometry (d) as a function of time from the start of observation for 
the 3.6 fim phase curve observations. Data have been binned into 
three minute intervals. In panels a and b solid horizontal lines indi- 
cate the pixel center while dashed horizontal lines indicate a pixel 
edge. Gaps in the data are due to spacecraft downlinks. Jumps in 
position, noise pixels, and relative flux after each downlink period 
are the result of stellar reacquisition. The pointing oscillations and 
long-term drift in the y direction are due to well-known instrumen- 
tal effects, and are present in all Spitzer observations. 

tions (Figures JT] and [2]), which is the result of intrapixel 
sensitivity variations aiscussed in the Appendix. 

2.2. 5.8 fim Photometry 

For our 5.8 /im data set, we determine the background 
level in each image using the same met hodo logy as pre- 
sented for the 3.6 and 4.5 /im data sets (S 2.1 ). The back- 
ground flux is 0.6% of the total flux in the science aper- 
ture. We tested both flux-weighted and Gaussian flts to 
the image methods of determining the stellar centroid. 
For the Gaussian fits to the image we first fit the image 
allowing both the x and y width of the Gaussian to vary. 
We then refit the same data set using a symmetric fixed 
width for the Gaussian that is the mean of the previous 
X and y widths. We flnd that Gaussian flts to the image 
give a lower standard deviation of the residuals in the 
final flts compared with flux- weighted centroids. This 
preference for the Gaussian centroiding method is likely 
the result of the short-time scale of these observations 
compared with the 3.6 and 4.5 /im observations and the 
less than 0.1 pixel change in the stellar centroid position 
over the course of the observation. 
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Figure 2. Measured x (a), y (b), noise pixels (c), and raw pho- 
tometry (d) as a function of time from the start of observation 
for the 4.5 fim phase curve observations. Data have been binned 
into three minute intervals. In panels a and b solid horizontal lines 
indicate the pixel center while dashed horizontal lines indicate a 
pixel edge; see Figure [l] for a complete description. 

We estimate the stellar flux from each background sub- 
tracted image using circular aperture photometry for a 
range of apertures sizes from 2.0 to 5.0 pixels in 0.25 
pixel increments. We find that an aperture size of 2.5 
pixels gives the lowest standard deviation of the residu- 
als in the final fits. We remove outliers in our data set 
more than 3. Oct away from a moving boxcar median 50 
points wide. Figure |3] show our resulting photometry for 
the 5.8 fxia observations. We test for possible intrapixel 
sensitivity variations using the methodology discussed in 
the Appendix, but find that including intrapixel sensi- 
tivity corrections actually degraded the precision of the 
final solution. 

2.3. 8.0 fj,m Photometry 

Unlike the 3.6, 4.5, and 5.8 data, the 8.0 /im data 
were observed in the full-array (256x256 pixels) mode of 
the IRAC instrument. Our photometry determines the 
flux from the star in a circular aperture of variable radius, 
as well as the background thermal emission surrounding 
the star. We calculate two values of the background, us- 
ing different spatial regions of the frame. First, we cal- 
culate the background that applies to the entire frame. 
Second, we isolate an annulus between 6 and 35 pixels 
from the star. In each region (entire frame, or annulus) 
we calculate the histogram of intensity values from the 
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Figure 3. Raw photometry as a function of time from the start 
of observation for the 5.8 fim secondary eclipse observation. Data 
have been binned into three minute intervals. Note the downward 
slope in the relative flux values as a function of time. 

pixels within the region, and we fit a Gaussian to that 
histogram. We use the centroid of that Gaussian as the 
adopted background value. This method has the advan- 
tage of being insensitive to 'hot' or otherwise discrepant 
values. We find that using the background values calcu- 
lated from the region near the star, as opposed to the 
entire frame, produces a lower scatter in the residuals in 
our final solution. 

We locate the center of the star by fitting a 2-D Gaus- 
sian to a 3x3-pixel median filtered version of the frame; 
we find that this method produces slightly better photo- 
metric precision than using a center-of-light algorithm. 
We estimate the stellar flux from each background sub- 
tracted image using circular aperture photometry for a 
range of apertures sizes from 2.0 to 5.0 pixels in 0.25 
pixel increments. We find that an aperture size of 3.5 
pixels gives the lowest standard deviation of the residu- 
als in the final fits. We remove outliers in our data set 
more than 3.0a away from a moving boxcar median 50 
points wide. Figure |4] show our resulting photometry for 
the 8.0 /im observations. We test for possible intrapixel 
sensitivity variations using the methodology discussed in 
the Appendix, but find that including intrapixel sensi- 
tivity corrections actually degraded the precision of the 
final solution. 

2.4. Flux Ramp Correction 

Our data exhibit a ramp-like increase, or decrease, in 
the flux values at the start of each observation and after 
each break in the 3.6 and 4.5 /im data for downlink. This 
ramp in the fiux values has bee n previously noted for 



IRAC 8 .0 /xm observation s (e.g. Knutson et al. (2007 



2009a); Agol et al. 



(|2010[)) as well as 3.6 and 4.5 iim 



observ ations (e.g. Beerer et al. (2011); Todorov et al. 



(2012 1). While the precise nature of this ramp is not 
well-understood, it can be at least partially attributed to 
thermal settling of the telescope at a new pointing, which 
contributes an additional drift in position. However, we 
see that the ramp often persists beyond this initial drift 
in position, and we therefore speculate that there is an 
additional component analogous to the effect observed 
at long wavelengths due to charge-trapping in the array. 
We tested several functional forms to describe this ramp. 
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Figure 4. Raw photometry as a function of time from tlie start 
of observation for the 8.0 fj-ra partial orbit phase curve observation. 
Data have been binned into three minute intervals. Note the strong 
exponential growth of the relative flux values as a function of time. 

including quadratic, logarithmic, and linear functions of 
time, but find that functional form that gives us the best 



correction is given by the formulation presented in Agol 
eFaLldMol): 



F'/F = 1 iflie^*/''^ iage 



-t/a4 



(1) 



where F' is the measured flux, F is the stellar flux in- 
cident on the array, t is the time from the start of the 
observation in days, and oi — are free parameters in 
the fit. 

We find that in the 3.6 and 4.5 ^m data the asymptotic 
shape of the ramp converges on timescales less than an 
hour for most of the data segments. Instead of including 
parameters in our fits for this ramp behavior, we instead 
elect to simply trim the first hour at the start of each 
observation and subsequent downlinks. This reduces the 
complexity of our fits while avoiding possible correlations 
between the shape of the underlying phase curve of HAT- 
P-2b and the ramp function. We find that the ramp per- 
sists beyond the one hour mark at the start of the 3.6 fim 
observations, which affects the shape of the first observed 
secondary eclipse and therefore include a ramp correction 
for that data segment. For the 5.8 /im data, we find that 
the fiux asymptotically decays from the start of the ob- 
servation (Figure pi). This decrease in the 5.8 /um fiux is 
also apparent in tne derived background flux values and 
has been pr eviously noted by studies such as 'Stevenson| 
et al. ( 201 1 ) . As seen in Figure |4) the 8.0 /im da ta exhibit 



the wefl know asymptotic increase in flux (e.g. Knutson 
et al.||2007[ [2009al [AjoTeFalllmol 



in order to select the best functional form for the ramp 
correction in each data set we use the Bayesian Informa- 
tion Criterion (BIC), defined as 



BIC ^ + klii{n) 



(2) 



where k is the degrees of freedom i n the fit and n is the 
total number of points in the fit (Liddle 2007). This 
allows us to determine if we are 'over-fitting' the data 
by including additional parameters to describe the ramp 
correction in our fits. We find that for the 3.6 /im data set 
using only a single exponential gives us the lowest BIC 
value. For the 8.0 /im data, using all the terms (ai — a^) 
in Equation[l]with no trimming of the data near the start 



of the observation gives us the lowest BIC value. For the 
5.8 /xm data, we find using a single decaying exponential 
gives the lowest BIC if we trim data within 30 minutes 
of the start of the observation. Trimming significantly 
more or less than than this amount ether gives a higher 
BIC or reduces the amount of out of eclipse data to less 
than the eclipse duration. 

2.5. Transit and Eclipse Fits 
We m odel our transit and eclip se events using the equa- 



tions of Mandel & Agol ( 2002 ) modified to account for 
the orbital eccentricity, e, and argument of periapse, u, 
of the IIAT-P-2 system. For an eccentric system the nor- 
malized separation of the planet and star centers, z, is 
given by 



r(t) 



y^l - (sini * sin (a; + f{t))y 



(3) 



where r{t) is the radial planet-star distance as a function 
of time, i?* is the stellar radius, i is the orbital inclination 
of the planet, and f(t) is the true anomaly as a function 
of time. The r{t)/Ri, term in Equation[3]is calculated as 



r(t) 



i?, l + ecos (/(<)) 



(4) 



The true anomaly angle (/) that appears in both Equa- 
tions[3|and |4lis determined using Kepler's equation (Mur- 
ray STDermott 1999|) and is a function of e, w, and the 
orbital period of the planet, P. Because of the degenera- 
cies that exist between e, lo and P in determining f{t) 
we elect to not use P as a free parame ter in our fi ts. In- 



stead we fix P to the value reported in Pal et al. ( |2010 ) 
5.6334729 d. We further minimize correlations in our so 
lutions for e and uj by solving for the Lagrangian orbital 
elements k = ecoso; and h = esincj. 

In addition to the parameters that define the orbit of 
HAT-P-2b, we solve for the fractional planeta ry radius, 
R^JR^. Because Rp/R^^ < 0.1 for HAT-P-2b ( |Pal et al. 
2010 Bakos et al. 2007b I, there exists a strong c orrela- 
tion between i a nd a/ Rj, in the transit solution (Winn 



e^e^2mh 



Pal 



2008|). We instead solve for the param- 

(|2007al 



eter s &^ and C/R* a s suggested by Bakos et al 
andlPal et al. (|2010|) where 



a 1 



R* 1 



cos I 



and 



R* 



a 2ti 



l + h 



R. p \/r^ vi - 



(5) 



(6) 



For the transit portion of the light curves we use four 
parameter nonlinear li mb-darkening coefficients for each 
bandpass calculated by Sing (2010), where we assume a 
stellar atmosphere with i e f f = 629 K, log{g) = 4.138, 
and [Fe/H]=-h0.14 (iPal et al.||2010[ ). For the secondary 
eclipse portion of the light curves we treat the planet 
as a uniform disk and scale the ingress a nd e gress to 
match the amplitude of the phase curve (j ]2.6[ ), which 
varies significantly over the duration of the eclipse and is 
therefore poorly approximated by a constant value. The 
secondary eclipse depth is defined by the average of the 
ingress and egress amplitudes. 
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For the 3.6 and 4.5 /im datasets, we use nine free pa- 
rameters to constrain the properties of the planetary or- 
bit, transit, and secondary eclipse. The 8.0 fxm obser- 
vations only include one secondary eclipse and therefore 
only require eight free parameters to constrain the same 
system properties. The 5.6 /im data set includes only the 
secondary eclipse portion of the light curve. We there- 
fore elect to only allow the secondary eclipse depth and 
timing to vary for the 5.6 /xm data and fix a/Ri,, i, e, uj, 
and Rp/Ri, to the average values from the 3.6, 4.5, and 
8.0 /zm data sets. 

2.6. Phase Curve Fits 

The functional form of the phase curve for a planet 
on an eccentric orbit is not well defined. Unlike close-in, 
tidally locked planets on circular orbits, eccentric planets 
experience time-variable heating and non-synchronous 
rotation ra tes. Previous studies b y [Langton fc Laughlin 
([2008 ) and [Cowan fc Agol[ ( [2011[ ) have investigated the- 
oretical light curves for plaiiets on eccentric orbits using 
two-dimensional hydrodynamic simulations and semi- 
analytic model atmospheres respectively. We also devel- 
oped a three-dimensional atmospheric model for HAT- 
P-2b that couples radiative and dynamical processes to 
further investigate possible phase curve functional forms 
that will be presented in a future paper. The functional 
forms for the phase curves described here all provide a 
reasonable fit to the theoret i cal light curves preseiit ed in 



Langton & Laughlin (2008L Cowan & Agol (20111, and 



trom our three-dimensional atmospheric model 

To first order, the flux from the planet is proportional 
to the inverse square of the distance between the planet 
and host star, r{t). This assumes that the planet has a 
constant albedo and responds instantaneously and uni- 
formly to changes in the incident stellar flux. We know 
that there must exist a lag in the peak of the incident stel- 
lar flux and the peak of the planet's tempe rature since 
atmospheric radiative timescales are flnite ([Langton fc 
LaughlinI 120081 liro fc Deming| [20T0l [Lewis ei al | [WHT 
Cowan & Agol 2011). Our first functional form for the 
phase variation of HAT-P-2b is a simple l/r(t)^ with a 
phase lag: 



= i^o + ci(l + cos(/-C2)) = 



(7) 



where / is the true anomaly and ci — C2 are free param- 
eters in the fit. The 1 -I- cos(/ — C2) is a simplified form 
of the denominator of Equation |4] for r{t). We also test 
a simpler form of Equation [7[ given by 



= Fo + ci cos(/-C2) 



(8) 



where ci — C2 are free parameters in the fit. This func- 
tional form of the phase curve is similar to a simple sine 
or cosine of the orbital phase angle, A or i^, used in the 
case of a circular orbit. 

We also find that a Lorentzian function of time pro- 
vides a reasonable representation of the expected shape 
of the orbital phase curve for HAT-P-2b. This is not 
surprising given that we expect the flux from the planet 
to vary as l/r(i)^ with a time lag between the mini- 
mum of the planetary distance and the peak of the plan- 
etary flux. We test both symmetric and asymmetric 
Lorentzian functions of the time from periapse passage. 



t, given by 



where 



Fit) = Fo 



Cl 



U(i)2 + 1 



U{t) ^ [t- C2)/C3 

in the symmetric case and 



[t - C2)/C4 



if t < C2; 

if t > C2 



(9) 
(10) 

(11) 



in the asymmetric case. In Equations[9[ |10[ and |ll| ci — C4 
are free parameters in the fit. The Lorentzian functional 
form for the phase curve is especially useful since the C2 
parameter gives the offset between the time of periapse 
passage and the peak of the planet's flux and the C3 — 
C4 parameters gives an estimate of relevant atmospheric 
timescales. 

We also find that the preferred functional form for the 
phase curve of planet s on circular orbits described in 
Cowan fc Agol (2008) provides a reasonable fit to our 
theoretical light curves if the orbital phase angle, ^, is 
replaced by the true anomaly, /. In this case 

F{e) = Fo + ci cos(6l) + C2 sin(e') (12) 
+C3 cos(26') -I- C4sin(26') 

where ci — C4 are free parameters in the fit and 6 = 
f + UJ + TT such that transit occurs a,t 9 = —Tr/2 and 
secondary eclipse occurs a.t 9 = n/2. 

In addition to the functional forms for the phase curve 
presented above, we also test a flat phase curve, F{t) = 
Fq, to make sure that we are not over- fitting the data 
with the phase curve parameters. In all cases we tie 
the Fq parameter to the secondary eclipse depth to give 
the appropriate average value of the phase curve during 
secondary eclipse. We also set any portion of the phase 
curve that falls below the secondary eclipse depth to zero. 
During secondary eclipse we no longer see flux from the 
planet, only the star, therefore the combined star and 
planetary flux should always be greater than or equal to 
the flux at the bottom of the secondary eclipse. 

We choose the optimal phase curve solution to be 
the one that gives the lowest BIC value (Equation [2| 
using a Levenbe rg-Marguardt non -linear least squares 
fit to the data ( [Markwardt] [2009[ ) . For the 3.6 and 
4.5 data sets, we also compare the BIC values 
from each of the three data segments separated by 
the downlink periods to make sure that the solution 
is robust at all orbital phases. For the 3.6, 4.5, and 
8.0 fim data sets, we fi nd t hat all of the functional 
forms presented Section 2.6 provide a significant im- 
provement in the BIC values over the 'flat' phase curve 
BIC values (BIC3 g ^m=l,282,575; BIC4 5 ^^=1,459, 141; 
BICg — 11,183). For the 3.6 /itm data we obtain the 
lowest BIC value with a functional form for the phase 
curve deflned by either Equation jSl or Equation [T2[ with 
C3 — C4 fixed at zero (BIC=1,278,521). This is not surpris- 
ing since basic trigonometric identities make Equation |12| 
equivalent to Equation [8[ if C3 — C4 can be assumed to be 
zero. For the 4.5 /xm aata we obtain the lowest BIC 
value with a fun ctional form for the phase curve given 
by Equation [12] with the ci and C3 terms fixed at zero 
(BIC=1,458,842). 
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The phase curve model given by Equation [12] also gives 
us a reasonable improvement in the BIG for our 8.0 fim 
data set. Ho wever, we find that the second harmonic in 
Equation [12] is highly degenerate with our ramp correc- 
tion at 8.0 /im such that we cannot reliably constrain the 
significance of this harmonic in the fit. We instead use 
the phase curve functional form given by Equation [7] 
which gives us the lowest BIG for the 8.0 /xm data set 
(BIG=10,945). For the 5.8 /im data, which only span 
ten hours, we find no statistical difference between solu- 
tions with and without parameterizations for planetary 
phase variations. 



2.7. Stellar Variability 



HAT-P-2 
(usini=20.7 

H & K lines of HAT-P-2 by [Knutson et al.[ (JMol do 



IS 

km 



a 



rapidly rotating F star 
Measurements of the Ga H 



not detect any significant emission in the line cores 
which would suggest a chromospherically quiet stellar 
host for HAT-P-2b. However, this indicator provides 
relatively weak constraints on chroniospheric activity for 
F stars, which have strong continuum emission in the 
wavelengths of the Ga II H & K lines. Ground-based 
monitoring of HAT-P-2 in the Stromgren b and y bands 
over a period of more than a year indicates that it 
varies by less than 0.13% at visible wavelengths. We 
would expect the star to vary by substantially less than 
this amount in the infrared, where the flux contrast 
of spots and other effects is correspondingly reduced. 
These observations rule out both periodic variability 
that could be associated with the rotation rate of 
HAT-P-2 or longer-term trends (G. Henry 2010, private 
communication) . Previous observations find the rotation 
rate of the star to be on the order of 3.8 d ( Winn et al. 
2007a ) based on the line-of-sight stellar rotation velocity 
(tJsinZi) and assuming sini,^ — 1, which is roughly 0.6 
times the orbital period of HAT-P-2b. We do not expect 
stellar variability to be significant on the timescales of 
our observations, but for the sake of completeness we 
also check this assumption directly using our Spitzer 
data. 

We employ two models for stellar variability. The first 
model is a simple linear function of time given by 



Foit) = dit 



(13) 



where t is measured from the predicted center of transit 
in each observation and di is a free parameter in the fit. 
The second model we test has the form 



Fo{t)=di sin((27r/d2)< - rfs) 



(14) 



where t is measured from the predicted center of transit 
in each observation and di — d^ are free parameters in the 
fit. Equation 14 attempts to capture stellar variability 
that is associated with star spots that rotate in and out 
of view with the d2 parameter representing the rotation 
rate of the star. 

We find that the linear model for stellar variability does 
not improve the of the fits significantly and in fact in- 
creases the BIG (Equatio n [2] ) . Inclusion of the stellar 
model given by Equation ]l4f does improve both the 
and BIG in our fits. However, we find that the solutions 
for the 'sine-curve' model for stellar variability are often 
degenerate with our models for the phase curve and the 



residual ramp at the start of the 3.6 /xm observations. We 
also find that the stellar rotation rate predicted by the 
d2 term in Equation 14 differs significantly between the 
3.6 and 4.5 /im observations with a rotation rate of 4.3 d 
preferred for the 3.6 /im data and 3.9 d preferred for the 
4.5 nm data. Although these predicted rotation periods 
are near to the expected 3.8 d rotation period of HAT- 
P-2, the amplitude of the predicted stellar variations in 
the mid-infrared seem spurious. The amplitude of the 
predicted stellar fiux variations at 3.6 and 4.5 /im are on 
the order of ~ 0.1%, which is comparable to our upper 
limit on stellar variability at visible wavelengths (0.13%). 
We would expect star spots to display much larger vari- 
ations in the visible than at mid-infrared wavelengths. 
The amplitude of these stellar variations are also similar 
to the amplitude of our predicted phase variations. We 
therefore conclude that our best-fit solutions for the stel- 
lar variability are physically implausible, and likely the 
result from a degeneracy with the other terms in our fits. 
As we have no convincing evidence for stellar variability, 
we assume that the star's flux is constant in our final fits. 

2.8. Radial Velocity Measurements 

Since the discovery of HAT-P-2b by |Bakos et al.| 
('2007b), the Galifornia Planet Search (GPS) team has 
continued to obtain regular radial velocity (RV) me asure 



ments of th is system using the HIRES instrument ( Vogt 
et al.[[1994|) on Keck. We used the GPS pipeline (see. 



e.g. ?) to measure precise RVs from the high-resolution 
spectra of HAT-P-2 using a superposed molecular iodine 
spectrum as a Doppler reference and point sprea d func- 
tion cahbrant. Previous studies of this system ( Bakos 



et al.l[2007bl [Winn et al.[l2007al [Pal et al.[[2010| ) have 
included HIRES radial velocity measurements through 



2007al 

ES radial velocity measurements through 
May 2008 (BJD = 2454603.932112). Here we present 
16 additional RV measurements from the HIRES instru- 
ment on Keck for a total of 71 data points spanning 6 
years (Tabl e[l]). We exclude me asurements obtained dur- 
ing transit (Winn "eral][2007a[ ) from our fits in order to 
avoid measurements attected by the Rossiter-McLaughlin 
effect. 

We simultaneously fit for the RV semi-amplitude (K), 
zero-point (7), and long term velocity drift (7) along with 
the orbital, transit, eclipse, and phase parameters in the 
3.6, 4.5, and 8.0 /im data sets separately. We also test for 
possible curvature in the RV signal (7), but find that our 
derived v alues f o r 7 w ere consistent with zero. As dis- 
cussed in Winn (2010), the transit to secondary eclipse 
timing strongly constrains the e cos uj term in our fits, but 
the e sin uj term is better constrained by the inclusion of 
these RV data. Rapidly rotating stars are known to have 
an increased scatter in their RV velocity distribut ion be- 
yond the reported internal error s (as discussed in [Bakos 



et al. 



a stel 



2007b[[Winn et al.[]2007a| . We therefore estimate 



ar Jitter term (crjitter) as described in Section [s] 



3. RESULTS 



We perform a simultaneous fit and calculate uncer- 
tainties for the relevant transit, secondary eclipse, phase 
curve, radial velocity, fiux ramp, and intrapixel sensitiv- 
ity correction parameters in our data set s using a M arkov 
Ghain Monte Garlo (MGMG) method ([Iw]p005[). For 



the 3.6, 4.5, and 8.0 iim data sets fit parameters include 
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Table 1 

Radial Velocities for HAT-P-2b from Keck 



BJD - 2400000 


RV 


Uncertainty 


(days) 


(m s^-"-) 


(m s^-*-) 


5o9ol. 777500 


-19.17 


7.46 


OvSyoz.oTlYUU 


-310.87 


7.64 


ooyoo.oi4ooo 


538.73 


T.oi 


50904.894980 


855.62 


7.94 


o4Uzd.D91ol9 


698.64 


7.94 




ana ^70 
696. /8 


7.72 


541o7. 104158 


684.36 


6.84 


o41o7. 159078 


717.17 


6.59 


04io5.UlDOOO 


757.70 


7.00 


o4ioo. io9DZz 


774.61 


6.80 


54189. 010o78 


651.86 


6.59 


o41o9.US8911 


635.16 


6.48 


cr^lQn it;'7'701 
04ioy.iO( (Zi 


619.70 


'7 on 


54ilD.959ci95 


722.59 


8.04 


54257. 7od4j1 


27.91 


5.79 


54Jo / . (OOO 1 1 


35.49 


D.iY 


'7(ir\'7r\n 
54io/. (bviOz 


24.33 


6.18 


54257.794116 


-11.70 


5.48 


54J5 / . (9D( /9 


-19.01 


5.02 


54io / . r994oz 


-20.94 


5.10 


cr ^ o 0001 /I o 

54257. oUzl45 


-23.48 


4.74 


54257. oU492d 


-21.89 


5.15 


54io / .oU mo4 


-35.65 


A m 
4.yi 


54257. olUo42 


-24.40 


5.08 


54257. olol4d 


-40.82 


5.22 


5425 / .olool r 


-37.08 


5.45 


5425 / .olo49U 


-23.98 


0.2i 


5425o.U2414d 


-313.00 


4.87 


54258. 0271d7 


-310.40 


4.89 


5425o.UoU292 


-311.56 


4.37 


54258. 0od278 


-326.04 


4.58 


5425o.U42doU 


-314.19 


5.10 


54258.045488 


1 0/^ 07 

-J^D.27 


4.92 


54258. 048ci9ci 


-342.10 


5.07 


54258. U514o^ 


-351.92 


4.88 


o42oo.Uo4ozo 


-356.80 


A V\A 

4.94 


54258. Uoolol 


-356.72 


A 

4.00 


cr ^ o cr o oic;i /i to 

54258. Uol472 


-360.49 


4.85 


54258. Uo4dod 


17/1 /I 7 

-J74.47 


4.73 


r: ,j o c o 0001 T O 

54258. U9911U 


-432.29 


5.27 


cr^ocro 1000 "D/^ 

54258. l(J2ooD 


-433.53 


4.97 


54258. 1(Jdd79 


-446.91 


4.74 


542 /y.o /Do9o 


371.51 


8.26 


O42o0.o2v5oo4 


137.15 


0.9i 


Cr-lOO/1 OTCTOO 

54294. S787U2 


TOO 70 

728.72 


6.61 


54oU4.8d49o2 


588.92 


0.07 


/I Qoc; 0*701 00 


737.55 


6.23 


54dOD.8D52lD 


731.39 


7.95 


54v>U7.912v579 


456.62 


6.39 


54v5oo.o12d19 


549. bU 


6.64 


cr-icr/ii:? 0001 7c 

54546.098175 


-684.07 


7.79 


cr^cr/17 iTcr 700 

54547. 1157UU 


536.15 


7.19 


54549. UoU4d8 


754.51 


'7 Of; 
7.20 


54dU2.91d5oU 


271.30 


a A'? 
0.47 


cr^^?oo oooiTO 

54dUo.9o2112 


662.48 


5.74 


548o9.1ddo95 


-369.59 


8.09 


rrcoiCC Q710Q1 

5oU1d.o i lUol 


a 70 


8.89 


crc^ocrr* o,icr/^oc 

55050. 945095 


-513.11 


8.45 


55405.742817 


606.08 


6.99 


55 ( O0.8 / 2995 


623.50 


7.70 


55704.836273 


377.02 


8.02 


55705.853346 


-535.77 


8.12 


55706.831645 


-257.19 


7.57 


55707.844353 


504.50 


7.68 


55808.759715 


323.97 


8.49 


55850.695496 


557.65 


7.37 


55932.161110 


-241.35 


9.48 


55945.128011 


595.28 


8.57 


55992.018498 


401.30 


10.08 


56147.753065 


553.56 


7.95 


56149.738596 


381.71 


8.37 



fo'^, C^jRi,, ecosw, esinw, Rp/Ri,, transit time, the sec- 
ondary eclipse depth(s), the phase function coefBcients 
ci — C4, a photometric noise term {(Jphot), and RV pa- 
rameters 7, 7, a jitter- For the 3.6 and 8.0 /im data 
sets we additionally fit for the ramp correction coeffi- 
cients fli — a4. The free parameters in the fit to the 
5.8 /im data set are the eclipse time, eclipse depth, ramp 
correction coefficients ai — 04, and aphot- The value of 
the stellar flux, Fq, is inherently accounted for with our 
intrapixel sensitivity correction method described in the 
Appendix. Because we do not apply intrapixel sensitivity 
corrections to our 5.8 and 8.0 /xm data sets we addition- 
ally fit for Fq in those cases. 

The only parameter held fixed in our analysis is the 
orbita l perio d (P), the value for which we take from Pal 



et al. (|2010|). All other orbital, planetary, phase, and 
ramp parameters are allowed to vary freely. Initial at- 
tempts at fitting just the 3.6 and 4.5 /im transits simul- 
taneously produce a value for P within la of the [Pal et al.| 
(20101 value, and we are therefore confident that adopt- 
ing the Pal et al. (2010) value for P has not introduced 
any signihcant errors into our analysis. 

We initially attempted to fit for the wavelength in- 
dependent parameters 6^, C/^*j ecosw, esincj, and 
simultaneously for the RV, 3.6, 4.5, and 8.0 /im data 
sets. However, given the size of our data sets (over 
2.5 million data points) and the time required to cre- 
ate our 'pixel-maps' for the 3.6 and 4.5 /im data (see 
Appendix), we found it computationally infeasible. It is 
possible that further improvements to our analysis code, 
including parallelization, could make the problem more 
computationally tractable. Such improvements are left 
for future iterations of our analysis methods. 

We plot the normalized time series for the 3.6, 4.5, and 
8.0 /im data sets after the best-fit intrapixel sensitivity 
variations and ramp corrections have been removed in 
Figure [51 The regions near secondary eclipse and transit 
for eacnbest fit solution are presented in Figures |6j [Tj [Sj 
and [9] for the 3.6, 4.5, 5.8, and 8.0 /im data sets respec- 
tively. We also present our best fit solution to the RV 
data in Figure [Tol 

We use a total of five independent chains with 10^ 
steps per chain in our MCMC analysis. Each chain is 
initialized at a position in parameter space determined 
by randomly perturbing the best-fit parameters from a 
Levenberg-Marquardt non-linear least squares fit to the 
data ( Markwardt 2009) . Instead of using a standard 
minimization scheme, we instead opt to maximize the log 
of the likelihood (L) given by 

log(L) = ^{- log(27rcr2) _ (data - model) V(2cr2)) , 

(15) 

where a is the relevant error term. This allows us to si- 
multaneously solve for the noise t erms gphot aud gutter 
with our other fit parameters (see Carter &: Winn||2009[ 
for further discussion of the maximum likelihood method 
as applied to exoplanet transit observations) . After each 
chain has reached 10^ steps, we find the the point where 
log(L) first surpasses the median log(L) value and dis- 
card all steps up to that point. We then combine the 
results from our independent chains and find the range 
about a median value that contains 68% of the values for 
a given parameter. We set our best-fit parameters equal 
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Table 2 

Global Fit Parameters 



Parameter 



3.6 ^im 



4.5 fim 



5.8 fj,m 



8.0 tim. 



Transit Parameters 

Rp/ 

\? 

in 

a/Rt 

Tc (BJD-2400000)'' 

C/R* (d-i) 

Tl4 (d)b 
Tl2 {d)b 



0.06821±0. 00075 



0.122 
87.37 



+0.066 
-0.078 
f-1.34 
-0.81 



9.53±0.35 
55288.84988±0.00060 
12.221±0.058 
0.1770±0.0011 
0.0128±0.0010 



o.o704i±o.ooo60 o.oegss'^ 

0.345±0.042 

84.91±0.47 85. 97^* 

8.28±0.24 8.70'^ 
55756. 42520±0. 00067 
12.286±0.057 
0.1813±0.0013 
0.0177±0.0012 



0.0668±0.0016 
0.2381°;?^? 

86.oj:i-o 

o oq+0.67 

54353.6911±0.0012 
12.21±0.11 
0.1789±0.0023 
0.0144±0.0021 



Secondary Eclipse Parameters 

Tl4 (d)!^ 
Ti2 {d)b 

I''* Eclipse Depth 
Tc (BJD-2400000)^ 
2'"' Eclipse Depth 
Tc (BJD-2400000)'' 



0.19177±0. 00029 

0.1550±0.0027 
0.01090±0. 00075 
0.080%±0.012% 
55284.2967±0.0014 

0.0993%±0.0090% 
55289. 9302±0.0014 



0.19310±0. 00025 
0.1651±0.0023 

0.01444±0. 00083 
0.1009%±0.0084% 
55751.8795±0.0011 
0.1057%±0.0090% 
55757.5130±0.0011 



n 071 +0-029% 

54906.8561+n nnl^ 



0.19253±0.00048 
0.1610±0.0043 
0.0121±0.0016 



0.1392%±0.0095% 
54354.7757±0.0022 



Orbital & RV Parameters 

e cos u) 

esiniii 

I n 

Tp (BJD-2400000) 

K (m s-l) 

7 (m s^-"-) 

7 (m s^-"- d^-'-) 



-0.50539±0. 00057 
-0.0741±0.0072 
0.51081±0. 00092 
188.34±0.80 
55289. 4734±0.0079 
927.0±5.9 
247.4±3.7 
-0.0890±0.0050 



-0.50301±0.00051 
-0.0729±0.0054 
0.50829±0. 00068 
188.25±0.62 
55757.05194± 0.0061 
923.0±5.8 
248.0±3.6 
-0.0881±0.0049 



0.50910=^ 
188. 09<* 



-0.50419±0.00088 
-0.0685± 0.0059 
0.50885±0.00097 
187.74±0.66 
54354.3109±0.0065 
923.1±6.0 
248.0±3.6 
-0.0886±0.0058 



Phase Curve Parameters 

Functional Form 

ci 

C2 

C3 

C4 

Amplitude 
Minimum Flux 
Minimum Flux Offset (h)'^ 
Maximum Flux 
Maximum Flux Offset (h)'^ 



Eq. (fL2l 
0.0379%±tTO015% 
0.0422%±0.0025% 

(fixed) 

(fixed) 
0.114%±0.010% 

"■""'Jl4/<'-o.oooi4% 

-18-361^0 
0.1139%±0.0089% 
4.39±0.28 



Eq. |fl2l 

(fijfea^ 
0.0293%±0.0042% 

(fixed) 
0.0163%± 0.0035% 
0.079%±0.013% 
n nq79%+o oo*6% 

0.0372 /0_Q gQgg5J 

6.71±0.43 

n llfi9(K+0.0089% 

5.84±0.39 



Flat 



Eq. 

0.0555%iftr.0032 
41.8°±2.7° 



0.1888%±0.0072% 
4.68±0.37 



Ramp Parameters 

Functional Form 

ai 

0,2 

03 

a4 



-0.00134l^gi°^ 
0.078l°:°|° 
(fixed) 
(fixed) 



None 



Eq. m 
+0.00683t^««o^^ 
0.095±0.015 



Eq. ITl 

-0 nn537+^™0088 

0195+°-°°^'^ 

-0.01765iroSS?? 
0.2812+0-0155 



Noise Parameters 

""phot 

CTjitter (m s-1) 



0.0042209±0.0000024 0.0057064± 0.0000032 0.015380±0. 000034 0.002164±0.000015 
26.0±2.1 25.7±2.2 ... 25.5±2.1 



^ We list all time in B,ID_UTC for consistency with other studies; to convert to BJD.TT add 66.184 s. 

^ Ti4 is the total transit or eclipse duration. T12 is the ingress duration, which equivalent to the egress duration (T34) to within 
error. 

^ Minimum ilux offset is measured relative to the center of transit time (Tc). Maximum flux offset is measured relative to the time 
of periapse passage (Tp). 

^ Represents average value from 3.6 /im, 4.5 fim and 8 fim analyses. Value held fixed in 5.8 /im analysis. 
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Figure 5. Final 3.6 (top), 4.5 fim (middle), and 8.0 fim photom- 
etry (filled circles) after correcting for intrapixel sensitivity varia- 
tions (3.6 and 4.5 /im) and the ramp-like behavior of the flux with 
time (3.6 and 8.0 /^m), binned into five minute intervals. The best- 
fit phase, transit and secondary eclipse curves are overplotted as a 
red line. The dashed line represents the stellar fiux level. 



to this median value and use this distribution to initially 
determine the la uncertainties in each of our parameters. 
For most parameters our error distribution was close to 
being symmetric about the median value. In the cases 
where the distribution was significantly asymmetric, we 
have noted both the positive and negative uncertainties 
in the parameter value. 

We find that there is ('red') correlated noise in our 
data even after the best fit intrapixel sensitivity varia- 
tions have been removed. To account for correlated noise 
in our data we first e mployed the wav e let-ba sed MCMC 
method described in Carter & Winn (2009). However, 
we found that our correlated noise could not be treated 
as a stationary noise parameter. Often we found that fits 
to the transit and secondary eclipse portions of our phase 
curves were degraded to introduce 'red' noise consistent 
with other portions of the light curve. As a result, we 
instead opt to use the 'residual-permutation' or 'prayer- 
bead' method to estimate the errors in our parame ters in 
the presence of correlated noise (see, for example Jenk 



ins et al..2002,,Southworth,2008,,Bean et al..2008, ,Winn 
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Figure 6. Best-fit transit (middle panel) and secondary eclipse 
(top and bottom panels) light curves (red lines) for the 3.6 ^im 
observations (black filled circles). The data have been binned by 
five minute intervals. Residuals for each of the fits are presented 
just below each transit or secondary eclipse event. The dashed red 
lines in the secondary eclipse panels show the best fit light curve 
corresponding to the other secondary eclipse. 



HAT-P-2b Phase Variations 



11 



1.0015 
^ 1.0010 



> 1.0005 



^ 1.0000 
0.9995 



0.0005 
0.0000 
-0.0005 



1.000 



> 0.998 
0.996 



0.0005 
0.0000 
-0.0005 



1.0015 
J 1.0010 

U- 

> 1.0005 



^ 1.0000 
0.9995 



0.0005 
0.0000 
-0.0005 



V' t I 

' . • . • . • 




— > i_ 



■ — X — ■ — T ■ M — X ■ — ■ — • — m—x — *T s a — ■ * = — ■ 



.* >• • •* s • 



-0.15 -0.10 -0.05 0.00 0.05 0.10 0.15 
Time From Eclipse/Transit Center (Earth Days) 



Figure 7. Best-fit transit (middle panel) and secondary eclipse 
(top and bottom panels) light curves (red lines) for the 4.5 ^m 
observations (black filled circles). The data have been binned by 
five minute intervals. Residuals for each of the fits are presented 
just below each transit or secondary eclipse event. The dashed red 
lines in the secondary eclipse panels show the best fit light curve 
corresponding to the other secondary eclipse. 
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Figure 8. Best-fit secondary eclipse light curve (red line) for the 
5.8 fim observations (black filled circles). The data have been 
binned by five minute intervals. Residuals to the fit are presented 
just below the secondary eclipse event. 
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Figure 9. Best-fit transit (top panel) and secondary eclipse (bot- 
tom panel) light curves (red lines) for the 8.0 lira observations 
(black filled circles). The data have been binned by five minute 
intervals. Residuals for each of the fits are presented just below 
each transit or secondary eclipse event. 
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Figure 10. Top Panel: RV measurements for HAT-P-2 presented 
in Table IT] (black circles) folded with an orbital period equal to 
5.633472g~^d along with the best fit RV solution (red line) with 
long term variation in the RV signal due to substellar object c 
removed. Open circles represent RV measurements taken after the 
completion of the analysis. Bottom Panel: Residuals between the 
RV measurements and best fit solution. Error bars on the data 
points include the best-fit stellar jitter term (crjittor)- 



et al. |2008). The 'prayer-bead' errors are typically 1.5- 
3x larger than the errors determined from our MCMC 
analysis. The largest increases in the uncertainty using 
the 'prayer-bead' method were for the planet-star ratio 
and secondary eclipse depths. In some cases the 'prayer- 
bead' errors are slightly smaller (~ 0.9x) than the errors 
from the MCMC analysis. In those cases we report the 
larger errors from the MCMC analysis. Table |2] presents 
the best-fit parameters for our data sets ana their la 
error bars. We find that our photometric errors, o-phot, 
are 1.05, 1.11, 1.11, and 1.15 times higher than the pre- 
dicted photon noise limit at 3.6, 4.5, 5.8, and 8.0 ^m 
respectively. 

4. DISCUSSION 

In the following sections we discuss the implications 
of these results for our understanding of the HAT-P- 
2 system and the atmospheric properties of IIAT-P-2b. 



HAT-P-2 system from 


Bakos et al. ( 


2007b), Winn et al. 


( 2007a 1, 


Loeillet et al 


( 


2008), ano 


I Pal et al. (2010|), 



velocity data. We also compare our results to predic- 
tions from one-dimensional radiative transfer and semi- 
analytic models of IIAT-P-2b's atmosphere. 

4.1. Orbital and RV Parameters 



We find that our estimates for the orbital and RV pa- 
rameters listed in Table [2] fall within the range of the val- 
ues from previous studies presented in Table |3] We note 
that there is often a more than 3a discrepancy between 
orbital parameters for the HAT-P-2 system presented in 
previous studies (Table [s]). We conclude that either pre- 
vious studies underestimate their error bars for the dis- 
crepant parameters, or the planet's orbital properties are 
varying in time. If we compare the orbital parameters we 
derived from the 3.6, 4.5, and 8.0 /im data sets, we find 
that our estimates are within 3a of each other. The pre- 
vious studies presented in Table [3] included only ground- 
based transit and radial velocity data. By the inclusion 
of secondary eclipse in our data we were able to improve 
the estimate of the orbital eccentricity of HAT-P-2 b by 
an order of magnitude over the value presented in |Pal| 
eraL] ( [2010| . 

From our orbital and RV parameters we can estimate 
the mass of HAT-P-2b using 



Mr. 



P" Gsini 



a 



Ri 



(16) 



where the orbital period (P) and stell ar radius [R*] ar e 



Pal et al.lpOlO ) 



assumed to be the values presented in 
5.6334729 ± 0.0000061 d and 1.64l°;°^i?0 respectively 
Values for the RV semi- amplitude {K), eccentricity (e), 
inclination (i), and normalized semi-major axis (a/Ri,) 
are taken as the error-weighted average of the values 
from the 3.6, 4.5, and 8.0 fim observations, which are 
presented in Table |4] We estimate the mass of HAT- 
P-2b to be 8.00 ± 0.97 Mj, which is within la of the 
previous estimates of Mp for HAT-P-2b (Table |3|. 

4.2. Linear RV Trend 

Our RV data spans a period of nearly six years, which 
allows us to test for long term trends in the RV signal. 
We find a non-zero value for the linear term in our RV 
fit (7), which indicates the presence of a second body 
orbiting in the system. If we assume that this second 
companion to HAT-P-2 ('c') is on a circular orbit and 
is much less massive than its host star, we can relate 
7 to the mass and orbi tal semi-major axi s of 'c' by the 
expression presented in Winn et al. ( 2009 ) 



7 ; 



GMc sin ic 



(17) 



Figure [TT] shows the range of values for sin ic and 
Qc that are allowed for the HAT-P-2 system given that 
7 = -0.0886 ± 0.0030. We know the orbital period of 'c' 
must be significantly longer than six years since we do 
not detect any significant curvature in our long term RV 
trend, so we set a lower limit of McSinic ~ 15 Mj and 
Qc ^ 10 AU based on an orbital period for 'c' of 24 years. 
We further employ adaptive optics imaging to search for 
'c' and set an upper limit on the values of McSinic and 

4.2.1. Adaptive Optics Imaging 

In an attempt to directly image the body responsible 
for causing the linear RV trend, we observed HAT-P- 
2 on May 29, 2012 using NIRC2 (PI: Keith Matthews) 
and the Keck II adaptive optics (AO) system at Mauna 
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Table 3 

Results from Previous HAT-P-2 Studies 



Parameter 



Bakos et al. ( 2007b i Winn et al. i 2007a I 



Loeillet et al. ( 2008 i 



Pal et al. I 2010 1 



Transit Parameters 

i n 

a/R* 

Tc (BJD-2400000) 
Ti4 (d) 
Ti2 (d) 



0.0684±0.0009 

>84.6 
q 77+1-10 

54212. 8563±0.0007 
0.177±0.002 
0.012±0.002 



0.0681±0.0036=- 
>86.8 
9.90±0.39=» 
54212.8565±0.0006 



U.0b89i_Q oQQgg 

",+0.85 
-0.93 



90.0"* 



10.28 



+0.12 
-0.19 



0.07227±0. 00061 

-1.12 
-0.87 



86.72^ 



8.99l«i? 
54387.49375±0.00074 
0.1787±0.0013 



0.0141 



+0.0015 
-0.0012 



Secondary Eclipse Parameters 
Tl4 (d) 

Tc (BJD-2400000) 



0.1847±0.0055'» 



0.1969±0.0040=' 



0.1896±0.0016'' 



0.1868±0.0019 
0.1650±0.0034 
54388. 546±0.011 



Orbital & RV Parameters 

P(d) 

e cos oj 

esinii) 

e 

CO n 

Tp (BJD-2400000) 
K (m s-i) 

Planetary Parameters 
Mp (Mj) 

Rp [Rj] 

Pp (g cm-^*) 
Qp (m s-2) 
a (AU) 



5.63341±0.00013 

0.520±0.010 
179.3±3.6 
54213.369±0.041 
1011±38 



9.04±0.50 
0.982l«;?g« 

7+44 
-16 



2271 
0.0677±0.0014 



5.63341 (fixed) 

0.501±0.007 
187.4±1.6 

883±57^ 

8.04±0.40 
0.98±0.04 
10.60±0.55'' 
207±20'' 
0.0681±0.0014^ 



5.63341 (fixed) 



0.5163 



+0.0025 
-0.0023 
-1.06 



189. 921^ 20 

(-0.0053 
-0.0030 

966.9±8.3 



54213.4798^ 



;.62 



+0.39 
0.55 
+0.039 



0-9511o.053 

2371^? 

7+0.0011 



0.0677^ 



0.0017 



5.6334729±0. 0000061 
-0.5152±0.0036 
-0.0441±0.0084 
0.5171±0.0033 
185.22±0.95 

983.9±17.2 



9.09±0.24 

1 ic:7+0.073 
l-l"^' -0.062 
7.29±1.12 

168±17 
0.06878±0.00068 



Noise Parameters 
o-jitter (m s~^) 



60 



31 



17 



Parameter value not directly quoted in reference, but calculated from quoted parameter values and errors 



Table 4 

HAT-P-2b parameters from a weighted 
average of the values from the 3.6, 4.5, and 
8.0 ^m fits 



Parameter 



Value 



Rp/ R* 

.'^transit 

(BJD) 

a/Ri, 

i n 



Esec 

E, 



(BJD) 

periapsc 

(BJD) 



Mp (Mj) 
Pp 



9p (m s" 
a (AU) 



0.06933±0. 00045 
2455288. 84923±0. 00037 
8.70±0.19 

0.50910±0. 00048 

188.09±0.39 
0.19253±0. 00018 
2455289. 93211±0. 00066 
2455289.4721±0.0038 
8.00±0.97 
1.106±0.061 
7.3±1.6 
162±27 
0.0663±0.0039 



Kea ( Wizinowich et ar||2000 ). Our observations consist 
of dithered images taicen with the K' (Ac — 2.12/ini) 
filter. Using the narrow camera setting, which pro- 
vides fine spatial samphng of the NIRC2 point-spread- 
function (10 mas pix~^), we acquired 9 frames each 
having 13.2 seconds of on-source integration time. The 
seeing was estimated to be 0.4" at visible wavelengths 
at the time of observations using AO wavefront sensor 
telemetry. We note that significant wind-shake degraded 
the quality of correction for some images but by only a 



marginal amount. 

The data were processed using standard techniques to 
correct for hot pixels, remove background radiation from 
the sky and instrument optics, flat-field t he a rray, and 
align and co-add individual frames. Figure [T2| shows the 
final reduced AO image along with the corresponding 
(IOct) contrast levels achieved as a function of angular 
separation. No companions were detected in raw or pro- 
cessed frames. 

We can use the limits from a non-detection to rule out 
the presence of companions as a function of Mc sin(ic) 
and Qc- Using HAT-P-2's parallax distance of 119 ± 8pc 
and estimat ed age of 2.6 ± 0.5 Gyr as determined by |Pal| 
et al. (20101 through a combined isochrone, light curve, 
and spectroscopic analysis, we find that our diffraction- 
limited observations are sensitive to companions on the 
hydrogen- fusing boundary for separations beyond ss 1". 
Interior to this region, the combination of imaging and 
RV data eliminates most low-mass stars, though late- 
type M-dwarf tertiaries located at ^ 40 AU could cause 
the long-term Doppler dri ft yet simultaneously evade di- 
rect detection (Figure 11) 



4.2.2. Orbital Evolution 

The likely presence of an M/L/T/Y dwarf at an orbital 
distance, Oc, of 10 to 40 AU from HAT-P-2b lends cre- 
dence to the possibility that HAT-P- 2b owes is c urrent 
orbit to a history of Kozai cycling (Kozai 1962) com- 
bined with long-running orbital decay generated by tidal 
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Figure 11. Top: RV variation as a function of BJD after sub- 
tracting tlie calculated variations due to HAT-P-2b. Red line shows 
our best-fit solution for the linear trend (7) in the data that re- 
sult from companion 'c'. Open circles represent RV measurements 
taken after the completion of the analysis, which conform with our 
measured linear trend. Bottom: Range of M sin i and semi-major 
axis (a) for companion 'c' (solid line) as estimated from the long 
term drift in the RV data (7). Dotted lines estimate the range in 
Msini and a given the error in 7. The shaded region gives the 
acceptable range of parameter space for companion 'c' given our 
upper and lower bounds on McS'mi vs. Oc (dashed lines). 



1998 



Wu & Murray 



2003 



Fab- 



friction ( Eggleton et al 

rycky fc Trcmaine 2007 ). At" first glance, this combina- 
tion ot Kozai Cycling and Tidal Friction (KCTF) seems 
more plausible than disk migration for producing the cur- 
rent orbital configuration of HAT-P-2b and indeed, long- 
distance disk migration for HAT-P-2b seems quite prob- 
lematic. The protostellar disk itself would have needed 
to be exceptionally massive and long lived, and an ad-hoc 
mechanism must be invoked to explain the large observed 
orbital eccentricity of HAT-P-2b. Furthermore, HAT- 
P-2b has planetary and orbital parameters that place 
it significantly outside the well-delineated population of 
'conventional' hot Jupiters with P ^ 3 days, M ~ Mjup, 
and e ~ 0. 

In the context of the KCTF process, HAT-P-2b is en- 
visioned to have formed well outside its current orbit, 
perhaps via gravitational instability in the original pro- 
tostellar disk. KCTF can operate if the mutual orbital in- 
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Figure 12. Left: Keck AO/NIRC2 K band image of the HAT- 
P-2 system. Image displayed using a square root scaling. Right: 
IOct contrast limit for companion detection as a function of angular 
separation from which define the upper bounds on Mc sini vs Oc. 

(or better, substantially larger) than the Kozai critical 
angle, ic = arccos [(3/5)1/2] _ (assuming > Mb, 
and an initially circular orbit for b). In the event that 
KCTF did operate, planet b experienced periodic cycling 
between successive states of high orbital eccentricity and 
high mutual inclination. The timescale for these cycles 
was 



2P2 M^-fMb + Mc 



37rPb 



Mr 



(1 



(18) 



With plausible values of Pc — 25 years, Cc — 0.5, and 
Mc — 20Afjup, TKozai — 300Kyr. By contrast, the gen- 
eral relativistic precession rate for HAT-P-2b is currently 



^GR — 



3G3/2(Af,-f Afb)3/2 



5/2 



(19) 



which generates a full 2tt circulation of the apsidal line 
in 20 Kyr (tgr)- Because TKozai ^ tgr, the magnitude 
of the Kozai cycling is strongly suppressed at present. 

To good approximation, tidal friction is currently the 
dominant contributor to the orbital evolution of HAT-P- 
2b, and it is plausible that the current state of the system 
has undergone dissipative evolution from an earlier epoch 
where rxozai = ^gr. If we assume that the tidal evolution 
has roughly conserved IIAT-P-2b's periastron distance, 
ab(l — eb) = 0.033AU, then in order for tko 



clination between companion 'c' and planet b was larger tgr implies that HAT-P-2b has evolved to its current 
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orbit from an orbit with an initial period, Po , given 
by 



Po„.„ = (1 - el) 



2 3 / 4 -fyW* / ^_G^ ^ 1 / 2 
C Mcdperi 



r', (20) 



or 



^o™>-.~33d(--^)( 



Pc ^^20Mjup^i/2 



25yr" 



(21) 



Additionally, the radial velocity-derived constraint on 

the unseen companion 'c' indicates that Mc ~ Pc^^, 
which allows us to simply write 



33d( 



Pc 

25 yr 



Nl/3 



(22) 



Given that direct imaging constrains the maximum pe- 
riod for companion 'c' to be of order 250 years, the KCTF 
scenario gives a lower limit on the possible initial periods 
for HAT-P-2b to be 30 days < Fo„„^ < 60 days. 

The KCTF process delivers planets into orbits in which 
the planetary orbital angular momentum and the stel- 
lar spin angular momentum are initially largely uncor- 
related. As mentioned earlier in the text, initial mea- 
surements of the projected spin-orbit misalignment an- 
gle, A, suggested that HAT-P-2b's orbit lies in the stellar 
equatorial plane . Rece ntly, however, an reassessment by 
Albrecht et~al] ( |2012[ ) reports A = 9° ± 5°, indicating 
a modest projected misalignment. In addition, HAT-P- 
2b's parent star is almost precisely on the Tog = 6250K 
boundary at which stars empirically appear to be abl e to 
maintain misalign ment over multi-Gyr time scales (Al- 
brecht et al.||2012|). 



4.3. Transit Timing Variations 

We calculate ephemeris for the HAT-P-2 system given 
the center of transit and eclipse times presented in Ta- 
ble [2] as: 

Tc{n) ^ T^{0) + n X P (23) 

where Tc is the predicted transit, eclipse, or periapse 
time and n is the number of orbits that have elapses 
since Tc(0). From our transit data we calculate Tc(0) = 
2455288.84923 ± 0.00037 BJD and P = 5.6334754 ± 
0.0000026 d, which is consisten t with the orbit al pe- 
riod for HAT-P-2b presented in iPal et al.l (|2010|). We 

It) 



calculate Tc{0) = 2455289.93211 ± 0.00066~BJb and 
P — 5.6334830 ± 0.0000086 d from our secondary eclipse 
data and Tc(0) = 2455289.4721 ± 0.0038 BJD and P = 
5.633479 ± 0.000026 d from our estimated times of peri- 
apse passage. The orbital periods estimated from the sec- 
ondary eclipse and periapse passage timings are within 
Icr of the value derived from our transit timings. Fig- 
ure [13] presents our transit, secondary eclipse, and peri- 
apse times compared to our derived constant ephemeris. 
We define the center of transit /eclipse to occur when 
the projected planet-star distance given by Equation [3] 
is minimized. The definition of transit /eclipse center is 
not always clearl y stated in studies for eccentric systems. 



Pal et al. (2010) estimate a difference between RV and 
photometric transit centers for the HAT-P-2 system of 
nearly two minutes. Some of the spread in the measured 
vs. predicted transit times could be accounted for by 
inconsistent definitions of transit center. 
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Figure 13. Observed minus calculated transit (top panel), peri- 
apse passage (middle), and secondary eclipse (bottom panel) times 
from data presented in Tables[2](squares) andp](triangles). Dashed 
lines indicate the la uncertainties in the predicted transit, periapse, 
and eclipse times. The variation that we see in the transit times 

The separate visits to HAT-P-2b for the 3.6, 4.5, 5.8, 
and 8.0 fj,m observations allow constraints to be placed on 
possible transit timing variations (TTV) for HAT-P-2b. 
Previous tran sit timing measurements for this system 



( ,Bakos et al 2007bl [Winn et al.||2007aj 



Pal et al. 

had relatively low timing resolution [At ^ 50 



20101 



and appeared to b e consistent with a constant ephemeris 



(Pal et al. 2010). Somewhat surprisingly, the Spitzer 
data trom orbits and 83 appear to be inconsistent with 
a constant ephemeris at the ^3.5cr level, with a typical 
deviation of 150 s. The deviations are anti-correlated 
between primary transit and secondary eclipse, and they 
switch sign on the two successive visits. This could be 
due to either an astrophysical cause, or an as-yet unmod- 
eled systematic error. In the sections below we consider 
two potential astrophysical explanations for the TTVs: 
a planetary satellite or an external perturber. 

4.3.1. TTVs generated by a planetary satellite 

IIAT-P-2b is expected to be in pseudo-synchronous ro- 
tation, in which its spin period is close to the orbital 
frequency in the vic inity of periapse. Given Porbit = 
5.63347d,|Hut|(|1981|)'s treatment gives 



5c:i 

16 



456* 



+ 



15e^ 



1 



■(1 



\3/2 /3g* 



-] 'Porbit = 1.89d. 

. . (24) 

In order to maintain orbital dynamical stability, a 
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prospective moon for HAT-P-2b must have have an or- 
bital semi-major axis, Osat, such that Ogat ^ FcritRwM, 
where i?Hiii is the Hill Sphere radius at periapse, 



i?Hm-api(l-e)(^)V3 



0.0042 AU, 



(25) 



and Fc„t < 0.5 (see, e.g. [Barnes fc 0'Brien| (|2002D for 



discussion of satellite stability tor planets with low orbital 
eccentricity). At distance Osat < FcritRwm from HAT-P- 
2b, the orbital period of the sateUite is Pgat = 0.385 d, 
which is substantially shorter than the Pspin — 1.89 d 
spin period of the planet. Tidal evolution will therefore 
cause the sat ellite's orbit to g radually spiral in to meet 
the planet ( |Sasaki et al.||20l"2 ). 

Based on the fairly uniform satcUitcs-to-planet mass 
ratios exhibited by the regular satellites of the Jovian 
planets in the solar system, and in the absence of any 
concrete information regarding exomoons, it is not unrea- 
sonable to expect Msat = Afpi/10^ 0.28M®. The char- 
acteristic time for a satellite's orbital decay, assuming an 
equilibrium theory of tides with a frequency-independ ent 
tidal quality factor, Qp, is (Murray fc Dermott]|1999 ) 



Qpi 



T- i3(^-»*^Hiii)''/'3^^^^^^^5^v g 



(^)i/^ (26) 



Adopting Qpi = 10^, and k2p = 0.1 we find r ~ 2 Gyr, 
which is comparable to the estimated stellar age, — 
2.7 Gyr ( |Pal et al.|20Tol ). The small value for the apsidal 
Love number, fe2p, stemming from the fact that HAT-P- 
2b at 8Mj up is quite cent rally condensed in comparison 
to J upiter Bodenheimer et al.| (see [2001 ) ; Batygin et al. 
(see 2009 tor a discussion of the relation between k2p and 
interior models) . Given the substantial range of possible 
values for Qpi, it would therefore not be unreasonable to 
find that HAT-P-2b harbors a fairly massive moon. 

A moon with the above qualities would produce transit 
timing variations of magnitude ( Kipping^^2009| ) 



1 /2 
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(27) 



where the geometric factor A appropriately amplifies or 
damps the transit timing variations in accordance with 
the eccentricity and orientation of the elliptical orbit of 
the planet about the parent star 



A = cos 



1 , — e cos(a;) , 
tan-i(- ^) 



- e sm uj 



2(1 -t- esinw) 
l-e2 



1 
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(28) 

Substituting the various values discussed above, we find 



^TTV^^t = 0.15 s. 



(29) 



which is several orders of magnitude too small to be of 
current interest. So while (somewhat surprisingly) it is 
distinctly possible that HAT-P-2b could currently harbor 
a massive satellite, any such satellite cannot be respon- 
sible for transit timing variations of the size that appear 
to be present. 

4.3.2. TTVs generated by an external perturber with a long 

period 



Perturbations from an as-yet undetected perturbing 
body present another potential explanation for the ap- 
parent transit timing variations. After the signature of 
HAT-P-2b's orbit has been removed from the existing 
RV observations, a secular acceleration that can be at- 
tributed to a distant companion remains. The constancy 
of the acceleration implies an orbital period for the com- 
panion of at least several years; for example, a body with 
mass, Mc ^ ISMjup, and semi-major axis, Oc = 10 AU, 
would suffice. For a configuration in which Oc S> api 
([Holman & Murray 2005), 



457r Mc 
T6~M^ 



(30) 



where ckg = api/(ac(l — ^c))- For a perturber with 



18Mj„ 



0.1s. 



0.3 and Ac = 10 AU, we find 
In addition, direct numerical in- 



M, 

tegrations indicate that for external perturbing planets 
that are consistent with the observed secular accelera- 
tion, the expected TTVs are invariably very small, and 
furthermore, do not exhibit the rapid variation shown by 
the timing reversal observed between orbit and orbit 
83. 

4.3.3. A low-mass resonant perturber? 

There are an effectively infinite number of stable orbits 
for perturbing bodies in mean-motion resonance with 
HAT-P-2b, and the diversity of such orbits is extended 
by HAT-P-2b's substantial eccentricity. A body in low- 
order resonance with HAT-P-2b can readily induce tran- 
sit timing variations of the magnitude that are appar- 
ently observed, and may be able to produce the curious 
structure exhibited by the timing variations of the sec- 
ondary eclipses and the primary transits. Exploratory 
calculations are currently underway to evaluate this pos- 
sibility. If we assume that the observed TTV are gen- 
erated by a gravitational perturber, this appears to be 
the most promising approach. However, we note that 
the significance of the reported transit timing variations 
is less than Aa. Further high-precision transit and/or 
secondary eclipse observations along with additional RV 
measurements would be needed to confirm the presence 
of these TTVs and constrain the properties of the per- 
turbing body. 

4.4. Transit and Secondary Eclipse Depths 

From our three transit observations we obtain planet- 
star radius ratios of 0.06821±0.00075, 0.07041±0.00060, 
and 0.0668±0.0016 at 3.6, 4.5, and 8.0 /Ltm respectively. 
Our estimates of the planet/star radius ratio, fl p/i?^, are 
significantly smaller than the value presented in Pal et al."] 
( 2010 ) , but a re fairly well ali g ned with th e values pre^ 
sented in Bakos et al. f2007b|, 
,Loeillet et al. (2008). Althoug 



Winn et a l. (2007a), and 



Pal et al. 



C2010 



mcor- 



2007b I 
ize the 



porate the 1 band observations from Bakos et ai 
they also include follow up observations that uti 
z and Stromgren b + y bandpasses. While the / and z 
bands probe a similar wavelength range (^ 0.8 — 0.9//m), 
the Stromgren b and y bands probe a slightly shorter 
wavelength range (~ 0.5/im) where atmospheric scatter- 
ing may become important. 

The difference between our values of Rp/Ri, at 3.6, 4.5, 
and 8.0 /im could point to enhanced atmospheric opac- 
ity in the 4.5 /zm bandpass due to CO, but our values 
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4 6 
Wavelength (|j.m) 

Figure 14. Secondary eclipse depths (squares) at 3.6, 4.5, 5.8, and 
8.0 /im compared to one- dimensio nal atmosphere models for HAT- 
P-2b's day side following [Burrows ct al. (2008) with the planet at 
a distance of 0.0478 A.U. i Ue models presented here incorporated 
values of the parameterized recirculation parameter (Pn) of both 
0.1 and 0.3 and allowed for the presence of a modest high-altitude 
absorber (k = 0.2 cm^ g"'^)- Filled circles represent the band inte- 
grated planet/star flux ratio for each of the IRAC bands. Models 
with an inversion best match data at 3.6, 4.5, and 8.0 fim. 



of Rp/Ri, differ by less than 2a. Given the large value 
of the average gravitational acceleration of HAT-P-2b 
(162±27 m s~^, Table W|, we would expect atmospheric 
scale heights to be smalX which supports our wavelength 
independent values for the planetary radius . If we as- 
sume a value of i?* — 



IMto fs^Q dPal et al. 



2010[ ), 

then we find an average radius value tor llAi-P-zb of 
Rp = 1.106 ± 0.061i?j. The radius of planet, given its 
mass, is in line with models of the thermal evolution of 
massive strongly irradiated planets. Planetary radii of 
1.13 to 1. 22 Rj for 8 Mj pla nets are expected for 1-10 
Gyr ages ( [Fortney et al]|2007| 



The estimates tor the secondary eclipse depths at 3.6, 
4.5, 5.8, and 8.0 fj,m presented here are the first sec- 
ondary eclipse measurements for HAT-P-2b. We com- 
pare these results to the predictions of atmospheric mod- 
els for this planet to better understand the atmospheric 
properties of HAT-P-2b. Figure [14] presents predicted 
day side emission as a function of wavelength from a one- 
dimensional radiative transfer model similar to those de- 



scribed in Burrows et al. (l2008|). These models assume 

pne 



a solar-metallicity atmosphere in local thermochemical 
equilibrium and include a parameterization for the re- 
distribution of heat from the day side to the night side 
of the planet (P„) and the possible presence of a high- 
altitude optical absorber (k) It is important to note that 
these one-dimensional models assume instantaneous ra- 
diative equilibrium and therefore do not account for the 
expected phase lag between the incident flux and plane- 
tary temperatures for eccentric planets. 

Models that also include an additional high-altitude 
optical absorber to produce a dayside inversion best 
match our 3.6, 4.5, and 8.0 /im data points. We find 
some difficulty in matching the 5.8 /im data point. We 
note that even with our best attempts to remove the 
ramp from this data that the secondary eclipse depth ap- 
pears to be slightly underestimated (Figure|8]). Given the 
large uncertainty in the 5.8 /im secondary eclipse depth 
it is within 2.5a of the ffux ratio predicted by our mod- 



els that include an inversion. Both the 4.5 and 8.0 //m 
secondary eclipse depths are more than lAa larger than 
those predicted by atmospheric models without an in- 
version. Even if the planet were assumed to be ^100 K 
warmer to account for the phase lag in planet-widc tem- 
peratures, models without an inversion would still under- 
estimate the planetary flux in the 4.5 and 8.0 /im bands. 

We calculate brightness temperatures in each band us- 
ing a PHOENIX stellar atmosphere model for the star 
(Hauschildt et al. 1999). The 3.6 /im eclipse depths 
have an average value ot 0.0996% ± 0.0072%, which cor- 
responds to a brightness temperature of 2392 ± 84 K 
We find an average 4.5 /im eclipse depth of 0.1031% ± 
0.0061% and a corresponding brightness temperature of 
2117 ± 65 K. Our eclipse depth at 5.8 /^m corresponds 
to a brightness temperature of 1613^^5q K , while our 
eclipse depth at 8.0 /im corresponds to a brightness tem- 
perature of 2258 ± 106 K. Our secondary eclipse mea- 
surements in the 3.6 and 4.5 /im bandpasses agree at 
the la level, and therefore do not provide any convinc- 
ing evidence for variability in the planet's ffux. Further 
observations of IIAT-P-2b at 3.6 and 4.5 /im could place 
better constraints on possible orbit-to-orbit variability in 
the planet's thermal structure. 

4.5. Phase Curve Fits 

The overall shape and timing of the maximum and 
minimum of the planetary ffux from our phase curves re- 
veals a great deal about the atmospheric properties of 
IIAT-P-2b. For planets on circular orbits, phase curve 
observations are generally related to day/night bright- 
ness contrasts on the planet. In the case of HAT-P-2b, 
the phase variations in the planetary flux are more in- 
dicative of thermal variations that result from the time- 
variable heating of the planet. In both the circular 
and eccentric cases the thermal phase amplitude and 
phase lag are determined by the radiative and advec- 
tive timescales of the planet. By comparing the 3.6, 4.5, 
and 8.0 fiui phase variations we can gain further insights 
into the thermal, wind, and possible chemical structure 
of HAT-P-2b's atmosphere. 

We ffnd that the peak of the observable planetary flux 
at 3.6 /im occurs 4.39 ±0.28 hours after periapse passage 
with a peak value of 0.1139% ± 0.0089%, which corre- 
sponds to a brightness temperature of 2556 ± 100 K. The 
exact timing and level of the minimum planetary flux at 
3.6 /im is much more uncertain. We flnd the the plan- 
etary flux at 3.6 /im drops below observable levels for a 
period of ^1 day 18.36^5 J* hours before transit. As such 
we can only put an upper limit of 1040 K on the mini- 
mum brightness temperature of IIAT-P-2b at 3.6 /im. 

The observed 4.5 /im HAT-P-2b phase curve exhibits 
a peak in the observable planetary flux 5.84 ±0.39 hours 
after periapse passage with a value of 0.1162%^° °°^^^^ 

and corresponding brightness temperature of 2255^^4 K. 
We detect emission from IIAT-P-2b over the entirety of 
its orbit at 4.5 /im with a minimum in the planet/star 
flux ratio of 0.0372%+°;°°^^^°, which corresponds to a 

brightness temperature of 1345^^33 K. This minimum in 
the planetary flux occurs 6.71 ± 0.43 hours after tran- 
sit. Roughly speaking, the shift of the minimum of the 
observed phase curve at 4.5 /im away from region be- 
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tween apoapse and transit points to a minimum in the 
planetary temperature that is shifted west from the anti- 
steUar longitude and/or that the nightside of the planet 
is still cooling even after the transit event. We would ex- 
pect the planet to have a temperature minimum east of 
the anti-stellar point near the sunrise terminator assum- 
ing a super-rot ating flow as shown for HD 189733b and 



HD 209458b in Showman et al. (2009). This dip in the 



planetary flux after transit is indeed puzzling and will be 
investigated further in a future study. 

Our 8.0 fim observations cover only the portion of 
HAT-P-2b's orbit between transit and secondary eclipse, 
so we can only constrain the behavior of the 8.0 fiui 
phase curve near the peak of the planetary flux. We 
find that the peak in the planetary flux at 8.0 /im oc- 
curs 4.64±0.33 hours after the predicted time for periapse 
passage. The maximum of the planetary flux at 8.0 /im 
is 0.1888 ± 0.0072%, which corresponds to a brightness 
temperature of 2806±79 K. 

It is significant that we obtain a good fit to the data 
using Equation 12 which is similar to the functional form 
used t o fit the hght curves of HD 1 89733b and WA SP- 
12b in Knutson et al. (2012) and Cowan et al. (2012a) re- 



spectively and produce a longitudmai map of the planet's 
thermal variations. The "longitudinal" direction of ther- 
mal phase maps, (j), must be understood to mean "local 
stellar zenith angle", Thermal phase variations probe 
the diurnal cycle, r(<I>). A tidally- locked planet has 1-to- 
1 correspondence between local stellar zenith angle and 
longitude (e.g.: = $), but phase maps can be made 



regard less of rotation rate (e.g., for Earth Cowan et al. 
2012b|. 



Nonetheless, there are two reasons that phase map- 
ping should not work for eccentric planets: 1) the time- 
variable incident flux makes the brightness map time- 
variable (T(<i>,i)), and 2) the time-variable orbital an- 
gular velocity of the planet dictates that longitudinally 
sinusoidal variations in brightness on the planet would 
not corre spon d to sinusoidal variations in the light curve. 
Equation 12 is sinusoidal in the true anomaly (/) rather 
than in time, which implicitly accounts for 2). 

The fact that Equation [T2] fits the phase variations of 
HAT-P-2b implies that thediurnal brightness profile of 
the planet is constant throughout the orbit: dT{^)/dt = 
0. This is entirely different from claiming that the longi- 
tudinal brightness profile of the planet is constant: since 
the planet is on an eccentric orbit, there is no fixed cor- 
respondence between longitude and stellar zenith angle. 
Rather, our data seem to indicate that the planet main- 
tains a constant brightness as a function of stellar zenith 
angle, in marked disagreement with expectations for such 
an eccentric planet. This is almost certainly due to a ge- 
ometric coincidence, however: HAT-P-2b shows us its 
day-side shortly after periapse. The day-side brightness 
of the planet likely changes throughout its orbit, but we 
are not privy to it. 

4.5.1. Interpreting the Flux Maximum 

We use the semi- analytic model developed in |Cowan 
Agol (2011) to interpret the peak amplitudes and phase 
lags of the planet's thermal brightness variations. This is 
essentially a two-dimensional, one-layer energy balance 
model where the user specifies Bond albedo, radiative 
timescale at the sub-stellar point at periapse, and the 



characteristic zonaQwind velocity. 

The Bond albedo is assumed to be 0.1 for all of the sim- 
ulations shown here. Measured albedos of h ot Jupiters 
range from 0.025 for Tr ES-2b ([Ki pping & S p]egel||20H' l 
to 0.32 for Kepler-7b ( jPemof^T et lirr[201 1] ) . Since 
HAT-P-2b is viewed equator-on, its albedo is degener- 
ate with poleward heat transport, which we do not ex- 
plicitly model. Increasing either albedo or meridionaj^ 
energy transport has the effect of uniformly decreasing 
the planet's flux throughout its orbit. 

The sub-stellar radiative timescale is spec ified at peri- 
apse, and is assumed to scale as Tr^ n oc (jShowman & 



Guillot|2002 Showman et al.|20T0 ), as one would expect 
for blackbody parcels of gas. The radiative relaxation 
time therefore varies throughout the orbit and is also a 
function of the location of a parcel of gas on the planet. 

Observationally, the zonal wind speeds on a gas giant 
like IIAT-P-2b are degenerate with its rotation rate. In 
our model we therefore specify the angular velocity of 
gas parcels in an inertial frame, rather t han in the usual 
rotating planet frame. By adopting the Hut ( 1981 1 pre- 
scription for pseudo-synchronous rotation ot binary stars, 
we convert the inertial angular frequency into an equa- 
torial zonal wind speed in the rotating planet's frame. If 
another prescripti on is more appropriate for th e planet's 
rotation rate (e.g. 
equatorial wind ve 
form offset. 



Ivanov fc Papaloizou 2007), then the 
ocities presented here are off by a uni- 



Both zonal wind speeds and albedo are assumed to 
be constant throughout the planet's orbit. The con- 
stant zonal wind velocity is likely the most problematic 
assumption given th at three-dimensional simulations of 



Kataria et al. 



([2012]) predict that zonal wind speeds at 
the mid-lK photosphere (pressures of 0.1-1 bar) change 
by tens of percent throughout the orbit of a hot Jupiter 
with an eccentricity of 0.5. It is also likely that the 
amount of equator-to-pole heat transport, which is de- 
generate with albedo in our model, will vary throughout 
HAT-P-2b's orbit. Our assumption of a constant albedo 
for the planet could also be limiting given that clouds 
could form near the apoapse of HAT-P-2b's or bit, then 
dissipate near periapse [Kane fc Gelino 2010). By fo- 
cusing on the region near periapse we limit the possible 
infiuence of temporal changes in albedo and wind speeds 
on our results. 

4.5.2. Model of HAT-P-2b and Circular Analog 

The top panels of Figure [T5[ shows how Tjad and Uwind 
affect the 4.5 peak flux ratio and phase lag for HAT- 
P-2b. The dependencies are qualitatively similar for the 
other two wavebands. For comparison, we also show 
in bottom panels of Figure [15] the same dependencies 
for a hypothetical circular twin to HAT-P-2b with semi- 
major axis fixed at the actual planets periapse separa- 
tion. There are a number of conclusions we can draw 
from these models: 

1. For circular planets, the radiative time scale and 
wind velocity are entirely degenerate: both the 
peak flux ratio and the phase lag depend solely 
on the pro duct Tradf wind Ci th e 'advective ef- 
ficiency' of Cowan & Agol (2011). Furthermore, 



^ zonal refers to the east /west direction 

^ meridional refers to the north/south direction 
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Figure 15. The dependence of peak flux ratio and phase lag at 4.5 fim on radiative timescale and zonal wind speed. The top panels are 
for HAT-P-2b; the bottom panels are for a hypothetical twin on a circular orbit with semi-major axis flxed at the periapse separation of 
HAT-P-2b. The plots were generated by computing 4.5 fim lightcurves on a 20 X 20 grid in r^ad a^nd '''wind using the energy balance model 
of,Cowan & Agol (201ll. 



the peak flux does not depend on the direction of 
zonal winds. The phase lag, on the other hand, 
would have the same amplitude but opposite sign 
for eastward vs. westward zonal winds. 

The peak flux ratio for the eccentric planet HAT- 
P-2b depends approximately on the advective ef- 
ficiency, e, as for circular planets, but the de- 
pendence is no longer monotonic. The peak flux 
ratio exhibits a maximum for radiative times of 
~ 6 hours and wind velocities of 2 km s~^. This 
is because eastward advection of heat brings the 
hot spot into view shortly afte r periapse, as dis- 
cussed in Cowan & Agol (20111. There is a local 
minimum in peak flux tor Tj-ad = 0, the radiative 
equilibrium case. The global minimum, however, 
still occurs in the limit of long radiative times and 
rapid winds, which results in a planet with zon- 
ally uniform, time-constant temperature, as in the 
circular case. 

Comparin g th e top-right and bottom-right panels 
of Figure 15 shows that zonal winds have quali- 
tatively the same effect, regardless of orbital ec- 
centricity: eastward winds make the peak flux oc- 
cur early, while westward winds cause a delay in 
the peak flux. The major difference between the 
two geometries is the phase-lag expected in the ab- 
sence of winds, the null hypothesis. For a circular. 



tidally-locked planet there is no phase lag in this 
limit, whereas the eccentric planet exhibits a large 
positive phase lag in the absence of winds. As a 
result of this different zero-point, the amplitude of 
phase lag from periapse decreases nearly monoton- 
ically for increasing e in the eccentric model, ex- 
actly the opposite behavior from a circular planet. 
Such counter-intuitive behavior is precisely why it 
is more useful to compare phase lags to the wind- 
less scenario ra ther than quote abso lute phase lags 
from periapse (Cowan fc Agol|[20lT|. 



4.5.3. Constraints on Model Parameters 

In Figure [16] we show exclusion diagrams in the Uwind 
vs. Trad plane for each waveband. Since we use a one-layer 
model, this can roughly be thought of as constraining the 
properties of the photospheres at each waveband, with 
the understanding that the mid-IR vertical contribution 
func tions span approximately a factor of ten in p ressure 
(e.g. |Knutson etal.|[2009bl [Showman et al.||2009[ ). Blue 
lines mTIgure show the 1 and '2a confidence mtervals 
from the peak flux value. Red lines show the same for the 
phase offset. Grayscale shows the combined confidence 
intervals. 

At 3.6 and 4.5 /im, the peak fluxes are consistent with 
a broad range of eastward zonal wind scenarios. Only 
very high values of e (top-right of plot) and most west- 
ward winds (bottom of plot) are excluded. At 8 /zm. 
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Figure 16. Exclusion diagrams for 3.6, 4.5, and 8.0 phase 
peaks. The blue lines show the one and two sigma confidence 
intervals based on the observed peak flux. Red lines show the 
same for the observed phase lag. Grayscale shows the combined 
constraint. The plots were generated by computing lightcurves on 
a 20 X 20 grid iii r^ari and fwind using the energy balance model of 
[Cowan fc Agol| | |2011[ |. 



only the highest peak flux values are favored, but even 
these are a very poor fit. If we adopt a Bond albedo of 
zero, the 8 /Ltm peak flux still disagrees with any mod- 
els by > 5(7. We therefore conclude that the high flux at 
8 /im is not due solely to atmospheric dynamics, but also 
to chemistry and the planet's vertical temperature pro- 
file. This waveband falls within a water vapor absorption 
feature and is therefore expected to originate higher up 
in the atmosphere than either the 3.6 and 4.5 /xm flux. 
The high flux in the 8.0 /im channel therefore suggests 
a temperature inversion, which is also supported by our 
measured secondary eclipse depths in Section [4. 4[ 

The constrai nts from the phase lag of peak flux (red 
lines in Figure 16) are stronger, albeit still degenerate. 



The data favor a narrow range of advective efficiencies 
('''radt^wind). The bottom panel of Figure |l6] gives the 
impression that the peak flux and phase lag at 8 /im 
combine to give a nice constraint on the model parame- 
ters, but the peak flux constraints should be taken with 
a grain of salt, as described above. It is important to 
note the decaying exponential-like trend of the preferred 
values of the zonal wind speeds with Trad in Figure |l6] 
that pairs short values for Tj-ad with higher values for 
^wind and longer values for Tj-ad with lower values for 
^wind. Three-dimensional models that consistently cou- 
ple radiative and advective processes will help to further 
constrain the relevant timescales in HAT-P-2b's atmo- 
sphere over the full range of its orbit. 

5. CONCLUSIONS 

In this paper we present the first secondary eclipse and 
phase curve observations for the eccentric hot Jupiter 
HAT-P-2b in the Spitzer 3.6, 4.5, 5.8, and 8.0 /im band- 
passes. These data include two full-orbit phase curves at 
3.6 and 4.5 /im, a partial-orbit phase curve at 8.0 /im, 
three transit events, and six secondary eclipse events. 
The timing between transit and secondary eclipse during 
a single planetary orbit combined with radial velocity 
measurements allows us to better constrain the eccen- 
tricity (e = 0.50910 ± 0.00048) and argument of periapse 
{uj = 188.09° ± 0.39°) of HAT-P-2b's orbit. Long-term 
linear trends in the RV data indicate the presence of a 
third body in the system. 

A comparison of our secondary eclipse depths with a 
one-dinicnsional model for the day-side emission of the 
planet suggests the presence of a day-side inversion in 
IIAT-P-2b's atmosphere. The timing and magnitude of 
the peak in the planetary flux at 3.6, 4.5, and 8.0 fiui 
are explained by a range of advective and radiative pa- 
rameters in our two-dimensional energy balance model, 
but suggest the presence of strong eastward equatorial 
winds (~ 2 — 6 km s~^) and short radiative timescales 
(~ 2 — 8 hours) at mid-IR photospheric levels near peri- 
apse. 

Further work is needed to fully explain our obser- 
vations of IIAT-P-2b. Three-dimensional atmospheric 
models that couple radiative and advective processes and 
include a range of compositions would help to further 
explain the phase variations we observe for IIAT-P-2b, 
especially outside of the orbital region near periapse. Ad- 
ditionally, low-resolution spectroscopy taken during sec- 
ondary eclipse would help to better constrain the atmo- 
spheric chemistry of IIAT-P-2b. Exoplanets on highly ec- 
centric orbits like HAT-P-2b present observers and mod- 



HAT-P-2b Phase Variations 



21 



elers with a unique opportunity to disentangle the contri- 
butions of radiative, advective, and chemical processes at 
work in hot Jupiter atmospheres. By refining our under- 
standing of exoplanets like HAT-P-2b we will be able to 
use that knowledge to better understand the atmospheric 
processes at work in other exoplanet atmospheres. 
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APPENDIX 
NOISE PIXEL PARAMETER 

The IRAC Point Response Function (PRE) results from the convolution of the stellar Point Spread Function (PSF) 
with the Detector Response Function (DRE) . The shape of the PRE is not constant and varies with the DRE at each 
stellar centroid position and with changes in the stellar PSF determined by the optics. The shape of the PREs differs 
substantially in each channel of the IRAC instrument as shown in Figure |17| Given this change in the shape of the 
PRE with wavelength, we expect that different methods to determine the stellar centroid position and correct for 
systematic effects may be needed in each of the four IRAC ban ds. Both the 3.6 a nd 4.5 /xm channels of the IRAC 
instrument are shortward of the Spitzer 5.5 /im diffraction limit ( Gehrz et al.|[2007 l and exhibit undersampled PREs 
whose shape changes as the stellar centroid moves from the center to the edge of a pixel causing intrapixel sensitivity 
variations. The PREs in the 5.8 and 8.0 /im channels of the IRAC instrument exhibit a more Gaussian-like shape that 
is better sampled by the detector resulting in only small intrapixel sensitivity variations. 

Ideally, we would like to account for these changes in the PRE in our photometric measurements. We cannot make a 
direct determination of the exact shape of the PRE as a function of pixel position, but we can measure the normalized 
effective-background area of the PRE (/3), which is also called noise pixels by the IRAC Instrument Team. If we 
assume the measured flux in a given pixel {Fi) is given by 

= FoP^ (Al) 
where Fq is the point source flux and Pi is PRE in the i*^ pixel. The noise pixel parameter, /3, is given by 

- (EF.r _ (EFoP.? _ iEP^? 



(A2) 



since Fq can be assumed to be constant. The numerator in Equation A2 is simply the square of the PRE volume (V) 
and the denominator i s the effective background a rea {(3). These quantities are related to the sharpness parameter, 
5*1, first introduced by MuUer & Buffington (1974) for constraining adaptive optics corrections by 



Si' 



(A3) 



The Si parameter describes the 'sharpness' of the PRE and can range from zero (flat stellar image) to one (all the 
stellar flux in the centr al pixel). 
As discussed in |Mighell ( 2005 ) , the Si parameter is related to the standard deviation (cr) of the PRE by 

1 



Si 



Cia^ 



(A4) 



where Ci is a constant that depends on the sampling of the PRE on the detector. From Equations |A4| and |A3| it can 
be shown that 



(A5) 



For both the 3.6 and 4.5 /xm data we measure /3 with the same circular aperture sizes used to determine the stellar 
centroid position. 

We find that /3 can serve two purposes in improving the signal-to-noise of the 3.6 ^m observations. First, using a 

circular aperture with a radius proportional to reduces the overall variations in raw unbinned flux from ~5% to 
~2%. Harris (1990) and Pritchet & Kline (1981) note that the optimal aperture radius for a circularly-symmetric 

Gaussian with a standard deviation of u is ro « 1.6 a, which is similar to our optimal aperture radius tq w ct ~ y/?. 
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4.5 [i\ 



8.0 |im 



Figure 17. Spitzer IRAC point response functions (PRFs) at 3.6, 4.5, 5.8, and 8.0 fim. The PRFs were generated by tiie IRAC team 
from bright calibrations stars observed at several epochs. The PRFs are displayed using a logarithmic scaling to highlight the differences 
in the PRFs in each bandpass. 



Second, we find that using y /? as an additional spatial constraint in the intrapixel sensitivity correction at 3.6 fim can 
improve the accuracy, as defined by the standard deviation of the residuals, in our final solution by ~l-2%. We find 
that using /3 as an additional constraint for the 4.5 /im observations does not significantly improve the accuracy of our 
results. Th is is not surprising since the IRAC 4.5 /im channel is closer to the Spiizer diffraction limit of 5.5 ^m (Gehrz 
et al. |2007i ). We also find that /3 does not vary with stellar centroid position in the 5.8 and 8.0 /im observations, which 
are longward of the Spitzer of the diffraction limit. 

fNTRAPfXEL SENSfTfVfTY CORRECTfON 

The 3.6 and 4.5 /im channels of the Spitzer IRAC instrument both exhibit variations in the measured flux that 
are strongly correlated with the position of the star on the detector (Figures |T] and [2]) . These flux variations are the 
result ofwell documented intrapixel sensitivity v ariations that are exacerbateob y an undersampled PRF (e.g.. Reach 



et al. (2005), Charbonneau et al. (2005 2008), Morales-Calderon et al. (|2006|), Knutson et al. (20081). The most 



common method used to correct mtrapixel sensitivity induced flux variations is to fit a polynomial function of the 



Morales-Calderon et al. 2006 Knutson 



stellar centroid position ( [Reach et al. 2005[ Charbonneau et al. 2005 2008[ ^ _^ 

et al.|[2008). This method works reasonably well on short time scaies (< 10 hours) wher e the variations in the s tellar 



centroid position are small (< 0.2 pixels). Recently, studies by Ballard et al. (2010) and Stevenson et al. 



implemented non-parametric corrections for intrapixel sensitivity variations by creating pixel 'maps 



w 



201 1| ) have 
iich give a 



smaller scatter in the residuals c ompared with param etric models in most cases. 

The pixel mapping method of Ballard et al. (2010) uses a Gaussian low-pass spatial filter to estimate the weighted 
sensitivity function given by 



where 



iVj -Vi? 



(Bl) 



(B2) 



i^Oj is the stellar flux measured in the jth frame, Xj/yj and Xi/yi are the stellar centroid position in the j'th and 
ith im age respectively, and Ux and Uy are the widths of the Gaussian fllter in the x and y directions. In the Ballard 
(2010) study, they assumed that all observations obtained outside of planetary transit (Jo,j) represented the 



et al. 



intrinsic stellar flux {Fq) convolved with the intrapixel sensitivity function. In our study, we must account for possible 
variations in the flux from t he II AT-P-2 system due to phase variations in the flux from IIAT-P-2b. This requires us 
to iteratively solve Equation Bl where Fq.j is determined at each step by removing a model for the planetary transit, 



eclipse, and phase variations trom the observ ed flux iF,, 



The main challenge in applying the Ballard et al. (2010) method to the HAT-P-2 data set is that as the number of 
data points, n, becomes large the time requ ired to compute Equation |B1| over the full data set becomes prohibitively 
long. We must iteratively solve Equation |B1| to constrain the phase variations of the planet, which requires the 
summation over n da ta points n times for each iteration. One solution is to bin the data into a manageable number of 
points as was done in Ballard et al. (2010). We find that binning the data degrades the precision of our final solution. 
Binning the data results in average values for the measured fiux and stellar centroid position that are not necessarily 
representative of the true va riati ons in the stellar flux due to intrapixel sensi tivity effects. We a lso flnd that the optimal 
Ux and Uy used in Equation 



B2 



varies with the stellar centroid position. ^Ballard et al. (2010) used fixed values for ax 
and Gy that empirically produces the lowest scatter. Using values for cr^ and Uy tJiat vary with centroid position gives 
us a lower scatter in the residuals compared to using fixed values for Ux and Uy. 



Here we present an enhanced version of the Ballard et al 
num ber of data points and optimized values of Ux and Uy 



wit 



B2 



(2010) pixel mapping method that allows for a large 
lout being computationally prohibitive. In Equation 
points that are outside ~ Q<Jx/y of the position of the ith data point will contribute negligibly to the weighted 
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Figure 18. Standard deviation of the residuals (circles) and intrapixel sen sitiv ity mapping time (triangles) as a function of the number 
of nearest neighbors included in the Gaussian weighting function (Equation |B4[ | for the 3.6 /im (top) and 4.5 /im (bottom) observations. 
The standard deviation of the residuals from our fits decreases rapidly as the number of nearest neighbors considered is increased from 5 
to 50 after wrhich the gains in the precision of the fit are negligible. The time required to compute the full pixel map increases more or less 
linearly with the number of nearest neighbors. We find that repeated iterations become too computationally expensive when more than 50 
nearest neighbors are considered. The computational time required per iteration is a multiple of the number of nearest neighbors and the 
1.2 million data points in each of the 3.6 and 4.5 ?m data sets. 



sensitivity function W{xi,yi). Given this fact, we reduce n in Equation Bl by only summing over a fixed number of 
nearest spatial neighbors. We determine the nearest neighbors to the ith hux measurement by the distance, r^, given 

by 



a{xj-Xi) +b{y.j~yi) +c 



3 1/2 



(B3) 



where a:,- 



and 



J are the position and noise pixel estimates for the jth. image. 

we ensure that the nearest neighbors to the zth flux measurement share the same 



B3 



By including /3 in Equation 

systematic effects. Although j3 is not strictly a spatial parameter, as discussed in the Appendix variations in j3 
incorporate systematics that affect the shape of the FRF including and beyond intrapixel sensitivity variations. Our 
use of /? is similar to the incorporation of pa rameters for the FSF width and elongation in the corr ectio n of systematic 
variations as discussed in Bakos et al. (2010) for HATnet data. The factors a, 6, and c in Equation B3 can be adjusted 
to give more weight to flux variations in a given spatial direction. For our 3.6 ^m. observations we find that the optimal 
value of yj^/b is 0.75 with a = c = 1. T his difference in th e a and b parameters accounts for the asymmetric shape of 
the IRAC FSF in the 3.6 /zm bandpass ( Gehrz et al.|2007 ). For our 4.5 /im observations we find that a = 6 = c= lis 
optimal although the results are similar if we assume c = 0. 

We calculate the Gaussian smoothing kernel K for the ith data point with respect to its jth nearest neighbor, Ki{j), 
according to 



(B4) 



and j, are determined by the standard deviation of the x, y, and ^^1"^ values for the n nearest 



where a^^i, (Jy,i, 

neighbors. This formulation for ax^i, 



ay^i, and cr^i/2 j gives a wider filter in poorly sampled regions and narrower filter 



in regions with higher spatial resolution. Since the x, y, and /3 vectors are independent of the fit parameters, we need 
only determine the n nearest neighbors and calculate Ki{n) once for use in our iterative fitting routines. 

When selecting the optimal number of nearest neighbors to use in this calculation, there is a direct tradeoff between 
increased precision and increased computational time. Figure [18] shows the change in the standard deviation of the 
residuals and the time required to compute Fqj for all j as a function of the number of nearest neighbors used. We find 
that using more than 50 nearest neighbors reduces the standard deviation of the residuals by less than 1%. We also 
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find a significant increase in tlie computational time required to use more tlian 20-50 ne ares t neiglibors. We tlierefore 
elect to limit the number of nearest neighbors considered in our calculation of Equation B4 to 50. 
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